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Selective  radio  frequency  (rf)  pulses  are  important 
in  both  NHR  spectroscopy  and  imaging.  Because  the  Bloch 
equation  that  dominates  the  dynamic  character  of  the 
magnetization  is  nonlinear,  the  design  of  the  rf  pulse 
shapes  capable  of  precisely  affecting  the  magnetization  in 
a  well-defined  frequency  range  is  not  a  simple  problem. 
Recently,  the  concept  of  optimal  control  has  been 
introduced  for  the  design  of  amplitude-modulated  rf  pulses 
as  a  time  function  to  provide  the  frequency  selectivity. 
Several  selective  pulse  shapes  have  been  developed  using 
these  methods;  however,  very  few  experimental  results 
that  are  consistent  with  the  numerical  calculations  have 
been  presented.  In  this  dissertation  we  will  give  a 


VII 


self-contained  and  detailed  derivation  and  explanation  of 
how  to  use  the  conjugate  gradient  method  to  solve  the 
optimal  control  problem.  We  will  present  the  simulation 
results  of  the  optimal  180  degree  and  90  degree  selective 
pulses  and  the  experimental  results  that  confirm  the 
numerical  calculations.  The  conclusion  is  that  the  optimal 
180  degree  and  90  degree  selective  rf  pulses  found  by  the 
conjugate  gradient  method  have  much  better  frequency 
selectivity  than  the  widely  accepted  sine  function  rf 
pulses.  Especially  interesting  we  will  prove  in 
experiments  and  simulations  that  the  responses  to  a  180 
degree  selective  pulse  will  be  the  same  no  matter  how  it 
is  used:  as  a  longitudinal  magnetization  inversion  rf 
pulse  or  as  a  transverse  magnetization  refocusing  rf  pulse. 
Thus,  the  conjugate  gradient  method  has  solved  the  problem 
of  designing  rf  pulse  shapes  for  frequency  selective 
excitation . 
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CHAPTER  I 
INTRODUCTION 

Selective  radio  frequency  (rf)  pulses  that  are  able 

to  excite  only  those  spins  lying  in  a  well-defined  slice 

are  an  essential  part  of  most  magnetic  resonance  imaging 

( MR  I )  rf  pulse  sequences.  Researchers  have  realized  for 

about  ten  years  that  a  conventional  sine  function  or  a 

sine  function  modified  by  a  Gaussian  function  rf  pulse  has 

severe  difficulty  exciting  only  those  spins  lying  in  a 

specific  slice  without  affecting  the  spins  outside  of  the 

slice,  especially  for  a  large  tip  angle.  Many  efforts  have 

been  made  to  improve  this  situation.  The  publications  in 

the  field  of  selective  rf  pulse  design  can  roughly  be 

divided  into  three  categories.  In  the  first  category 

an  analytical  method  is  used  to  analyze  and  to  design  rf 

pulses.  The  representative  authors  are  G.  A.  Morris  and 

R.  Freeman  (1),  D.  I.  Hoult  (2)  and  W.  S.  Warren  (3). 

Morris  and  Freeman  treated  the  Bloch  equation  as  a  linear 

differential  equation  and  applied  the  convolution  theorem 

to  find  the  selective  pulse  sequence  called  DANTE  that 

consists  of  a  series  of  non-selective  rf  pulses.  Since  the 

Bloch  equation  is  essentially  not  linear  in  rf  pulses  and 

DANTE  does  not  combine,  a  magnetic  field  gradient  that  can 

be  used  to  discriminate  the  spatial  difference,  the  pulse 


sequence,  DANTE,  has  few  applications  in  MRI.  Hoult 
applied  perturbation  theory  to  the  Bloch  equation  and 
analyzed  the  rf  pulses.  Warren  applied  perturbation  theory 
to  the  liarniltonian  of  the  magnetic  resonance  spin  system. 
Since  selective  90  degree  and  ISO  degree  rf  pulses  cannot 
be  treated  as  perturbations,  the  results  of  Hoult  and 
Warren  must  be  very  rough  and  restricted  in  applications. 
In  fact,  the  selective  rf  pulses  that  they  found  have  a 
worse  slice  selectivity  than  that  of  the  optimal  selective 
rf  pulses  found  by  the  conjugate  gradient  method  that  will 
be  discussed  in  the  following  chapters. 

In  the  second  category  an  analytical  method  combined 
with  intuitions  is  used  to  analyze  and  to  design  the  rf 
pulses.  The  representative  authors  are  D.  G.  Nishimura  (4) 
and  H.  Yan  and  J.  C.  Gore  (5).  They  used  two  or  more 
amplitude-modulated  rf  pulses  each  of  which  combined  a 
specific  magnetic  field  gradient  and  constructed  a  single 
selective  rf  pulse.  The  selectivity  of  their  pulses  is 
worse  than  that  of  the  optimal  rf  pulses  and  the 
complication  in  practice  to  implement  their  pulses  is 
higher  than  the  single  optimal  rf  pulses.  Because  of  these 
reasons  Nishimura  never  uses  his  own  pulses  in  his  later 
publications  and  instead  he  uses  the  optimal  rf  pulses. 

In  the  third  category  a  synthetic  method  is  used  to 
design  the  selective  rf  pulse.  "Synthetic"  means  here  that 
one  can  start  from  a  desired  slice  profile  and  try  to  find 
an  actual  slice  profile  that  will  approach  the  desired 


slice  profile  as  closely  as  possible  by  changing  the  rf 
pulse  shape  or  phase.  The  representative  authors  are 
M.  S.  Silver  et  al.  (6),  A.  1! .  Lent  and  M.  R.  Kritzer  (7), 
S.  Conolly  et  al.  (8)  and  M.  O'Donnell  and  W.  J.  Adams  (0). 
Silver  et  al.  proposed  the  use  of  a  complex  hyperbolic 
secant  excitation  function  as  a  180  degree  inversion  rf 
pulse.  Their  result  performs  well  for  magnetization 
inversion,  but  they  have  not  discussed  the  application  of 
their  pulses  in  the  transverse  magnetization  refocusing 
and  the  pulse  length  is  about  25  ms  which  is  too  long  for 
MRI.  Thus,  their  pulses  are  not  suitable  for  MRI.  Lent  and 
Kritzer,  Conolly  et  al.  and  O'Donnell  and  Adams  found  that 
the  design  of  selective  rf  pulses  can  be  reduced  to  an 
optimal  control  problem,  or  a  dynamic  optimization  problem. 
Researchers  using  different  methods  to  solve  the  optimal 
control  problem  found  different  results  and  very  few 
experimental  results  that  are  consistent  with  the 
calculations  have  been  published. 

The  article  of  Conolly  et  al.  is  especially 
interesting.  We  try  to  give  a  little  more  detailed 
comments  about  their  article.  Conolly  et  al.  had  two  goals 
in  their  article.  First  is  to  give  an  efficient  algorithm 
to  find  the  optimal  selective  rf  pulses.  The  method  used 
by  Conolly  et  al.  to  optimize  the  objective  function  is 
the  steepest  descent  method.  The  steepest  descent  method 
is  a  less  efficient  convergent  iteration  procedure  than 
the  conjugate  gradient  method  (10).  The  results  found  by 


Conolly  et  al.  are  not  quite  correct  because  they  did  not 
use  their  optimal  rf  pulses  in  their  recent  article 
discussing  the  MR  angiography  by  selective  inversion 
recovery  (11)  and  instead  they  use  a  optimal  pulse  that  is 
quite  similar  with  our  optimal  pulse  (12).  A  reason  that 
Conolly  et  al.  could  not  find  the  correct  results  may  be 
the  slow  convergent  rate  of  the  steepest  descent  method. 
Because  of  the  slow  convergent  rate  of  the  algorithm  the 
authors  had  no  good  chance  to  conduct  many  calculations 
and  to  find  the  correct  parameters  leading  to  the  correct 
results  by  using  the  trial  and  error  method.  Another 
reason  that  Connolly  et  al.  could  not  find  the  correct 
results  may  be  that  they  chose  an  incorrect  desired  slice 
profile  because  they  believed  that  any  desired  slice 
profile  could  be  reached  by  the  optimization  procedure  and 
they  did  not  realize  the  relation  between  the  optimal  rf 
pulse  shape  and  the  desired  slice  profile. 

The  second  goal  of  Conolly  et  al.  was  to  prove  the 
existence  of  practical  selective  90  degree  and  180  degree 
amplitude-modulated  rf  pulses  that  can  cause  the  actual 
magnetization  to  have  any  desired  magnetization 
distribution  by  proving  the  controllability  of  the 
magnetic  resonance  imaging  spin  system.  Here,  "practical" 
means  that  the  pulse  length  cannot  exceed  the  maximum  rf 
pulse  length  used  in  actual  magnetic  resonance  imaging  and 
the  pulse  amplitude  cannot  exceed  that  which  an  actual  rf 
pulse  amplifier  can  supply.  Thus,  the  control  problem  is 


one  with  a  control  function,  which  is  the  rf  pulse  in  MR  I, 
constrained  in  both  pulse  amplitude  and  pulse  length.  The 
authors  did  prove  that  the  magnetic  resonance  imaging  spin 
system  is  controllable,  however,  the  authors  did  not 
mention  that  transferring  to  any  desired  distribution  may 
require  a  no n- practical  pulse  amplitude  corresponding  to  a 
practical  pulse  length  or  a  non-practical  pulse  length 
corresponding  to  a  practical  pulse  amplitude,  even  though 
that  length  is  finite  and  definite.  In  other  words,  that 
any  desired  distribution  is  reachable  was  not  proven  for 
the  control  function  constrained  in  both  amplitude  and 
length . 

Conolly  et  a  1 .  realized  that  the  pulse  length  may  be 
sufficiently  long  for  the  proof  of  the  controllability, 
but  the  authors  did  not  point  out  that  a  sufficiently 
long  pulse  length  may  be  not  practical.  In  the  discussion 
of  minimum  distance,  the  authors  realized  that  we  could 
trade  off  some  slice  definition  for  weaker  pulses,  but  the 
authors  did  not  point  out  that  a  special  slice  definition 
may  lead  to  a  stronger  pulse  and  the  amplitude  of  the 
pulse  may  be  not  practical.  The  practical  and  reasonable 
pulse  length  and  pulse  amplitude  have  significant 
importance  for  practical  rf  pulse  design.  Failing  to 
consider  these  restrictions  in  the  proof  proves  nothing 
for  the  goal  above.  The  controllability  of  the  spin  system 
without  considering  these  restrictions  is  only  of  academic 
interest.  These  restrictions  are  not  clearly  indicated  by 


the  paper  and  the  conclusion  of  the  controllability  is 
misleading  for  practical  rf  pulse  design.  Especially,  the 
reader  may  be  misled  into  thinking  that  any  desired 
distribution  is  reachable  even  with  the  contrained  pulse 
length  and  pulse  amplitude  because  the  authors  had 
discussed  the  finite  pulse  length  and  the  pulse  magnitude 
limit  in  several  places  in  the  paper.  Actually,  the 
conclusion  of  the  paper  that  the  optimal  control  design 
approach  allows  for  inversion  of  the  Bloch  equation  for 
any  desired  distribution  is  obviously  not  correct  for  a 
practical  rf  pulse.  For  example,  one  cannot  force  a  slice 
edge  which  corresponds  to  a  certain  slice  thickness  to 
have  every  finite  slope  by  using  a  single  rf  pulse  of 
practical  pulse  length  and  practical  pulse  amplitude.  One 
can  only  say  that  in  practice  some  desired  magnetization 
distributions  can  be  reached  but  others  can  not. 

We  should  also  point  out  that  in  the  article  of 
Conolly  et  al.  the  conclusion  that  all  optimal  90  degree 
pulses  have  the  same  shape  is  not  correct.  This  mistake  is 
because  the  authors  did  not  pay  any  attention  to  the 
relation  between  the  optimal  rf  pulse  shape  and  the 
profile  requirement  for  a  desired  slice,  for  example,  the 
requirement  for  a  definite  slice  edge  slope  of  a  desired 
slice  profile  with  a  certain  slice  thickness.  A  different 
requirement  for  the  slice  edge  slope  leads  to  a  different 
optimal  rf  pulse  shape.  This  conclusion  can  also  be 
applied  to  the  optimal  180  degree  rf  pulses.  Ke  will  give 


a  detailed  discussion  of  this  problem  in  the  Chapter  IV 
of  this  dissertation. 

We  should  mention  that  the  reason  why  an  actual 
slice  profile  obtained  by  an  optimal  rf  pulse  always  has 
some  errors  when  comparing  with  the  desired  slice  profile 
is  because  any  desired  distribution  is  not  reachable  by  a 
practical  rf  pulse,  not  because  of  the  formulation  of 
minimum  distance  as  claimed  by  the  Conolly  et  a 1 .  in  their 
article.  In  practical  rf  pulse  design  the  squared  error 
objective  function  cannot  achieve  zero  by  simply 
increasing  the  pulse  length  of  a  practical  rf  pulse. 
Thus,  the  main  results  of  the  article  of  Conolly  et  a 1 . 
are  basically  incorrect  even  though  Connolly  et  a  1 .  are 
one  of  the  first  to  use  the  optimal  control  concept  to 
design  the  selective  rf  pulses. 

In  this  dissertation,  we  try  to  give  a  detailed  and 
self-contained  derivation  and  explanation  of  the 
application  of  the  conjugate  gradient  method  (10,12,13)  to 
design  the  selective  rf  pulses.  We  will  discuss  the 
experimental  procedure  and  give  the  experimental  results 
that  are  consistent  with  our  theoretical  calculations. 
Comparing  the  experimental  results  with  the  calculated 
results,  one  can  conclude  that  the  optimal  control  method 
has  solved  the  selective  rf  pulse  design  problem  that  has 
existed  since  the  appearance  of  the  magnetic  resonance 
imaging.  The  results  presented  in  this  dissertation  can  be 
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considered  as  an  expansion  of  our  previous  results 
(12,14,15). 

Chapter  II  will  be  devoted  to  show  how  the  design  of 
a  180  degree  selective  rf  pulse  can  be  reduced  to  an 
optimal  control  problem.  Chapter  III  will  give  a  detailed 
derivation  of  how  to  use  the  conjugate  gradient  method  to 
solve  the  optimal  problen.  Chapter  IV  will  give  the 
calculated  results  for  several  special  cases  of  180 
degree  selective  rf  pulses  and  compare  the  calculations 
with  the  experiments.  Chapter  V  will  give  a  detailed 
derivation  and  procedure  of  how  to  use  the  conjugate 
gradient  method  to  design  the  90  degree  selective  rf 
pulses.  Chapter  VI  will  discuss  several  special  cases  of 
the  90  degree  selective  rf  pulse  and  compare  the 
calculations  with  the  experiments. 

In  the  discussion  in  the  folio w in g  chapters,  some 
variables  may  be  ignored  for  brevity  from  their  relevant 
functions  if  they  are  evident  from  the  context.  For 
example,  the  magnetization  component  Mz(t,l,U)  has  three 
variables,  t,  1,  and  U,  and  sometimes  we  just  write 
i i z ( t , 1 , U )  as  M z ( t )  or  M z ( t , 1 ) ,  where  t  is  time,  1  is  an 
integer  for  frequency  step  numbers  and  U  is  the  rf  pulse 
used  as  an  input  to  the  magnetic  resonance  spin  system. 


CHAPTER  II 
DESIGNING  180  DEGREE  SELECTIVE  INVERSION  RF  PULSES 
IS  AN  OPTIMAL  CONTROL  PROBLEM 

The  dynamic  equation  of  the  spin  magnetization  is  the 

Dloch  equation 


dM      _   _    Ml*i  +  M2*j     (M3  -  M0)*k 

=  r*yt   x  II 

dt  T2  Tl 


[1] 


where  Ml,  M2  and  M2  are  x,  y  and  z  components  of  the  spin 
magnetization  vector  M ;  MO  is  the  equilibrium  value  of  the 
magnetization;  K  is  the  applied  magnetic  field;  Tl  and  T2 
are  the  longitudinal  and  transverse  relaxation  times, 
respectively;  i,  j  and  k  are  the  unit  vectors  of  the  x,  y 
and  z  axes;  the  asterisk,  * ,  means  multiplication;  the 
bar,  -,  means  vector;  r  is  the  gyro magnetic  ratio. 

Because  the  Bloch  equation,  Eq.  [ 1  ]  ,  is  very 
complicated  for  designing  a  selective  rf  pulse,  Me    only 
consider  the  situations  in  which  the  relaxation  and 
diffusion  effects  can  be  ignored  when  the  selective  rf 
pulse  is  acting  in  a  pulse  sequence  and  Eq.  [1]  can  be 
reduced  to  the  following  equation: 


dM 

=  r*M  x  II. 

dt 


[2] 


The  applied  magnetic  field  can  be  expressed  as 

H  =  HO  +  Gz*z*k  +  Hl(t)*I  [3] 

where  HO  is  the  main  magnetic  field  which  is  assumed  in 
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N  X  I 


I'  = 


the  z  direction;  Gz  is  the  magnetic  field  gradient  applied 
along  the  z  direction;  Ill(t)  is  the  selective  rf  pulse 
applied  along  the  x  direction  on  resonance.  If  one 
consider  the  problem  in  the  rotating  frame,  the  applied 
magnetic  field  in  the  rotating  frame  is 

H  =  Gz*z*k  +  Hl(t)*I.  [4] 

Because  of  the  folio  wing  relation: 

i     j     k 
Ml     M2     M3 
HI     0    Gz*z 
=  M2*Gz*z*i  -  Ml*Gz*z*"J  +  H1*M3*J  +  Hl*M2*k, 

[5] 
the  Bloch  equation  without  relaxation  in  the  rotating 
frame  is 


[6] 


where  H(z)  =  z  *  G  z  is  the  magnetic  field  that  is  caused  by 
the  magnetic  field  gradient  Gz;  U(t)  =  -r*Hl(t)  ,  Ill(t)  is 
the  selective  rf  pulse  and  it  can  be  considered  as  an 
input  to  the  Bloch  equation.  For  convenience  of  discussion, 
we  can  call  either  Hl(t)  or  U(t)  as  the  selective  rf  pulse. 
For  convenience  of  the  following  discussion,  one  can  write 
Ml,  M2  and  M3  more  detailed  as  Ml(t,l,U),  M2(t,l,U)  and 
M3(t,l,U),  respectively. 

An  ideal  selective  rf  pulse  will  excite  only  those 
spins  lying  in  a  well-defined  slice  without  any  effect  to 


'   111 

0 

r  *  H  (  z  ) 

0 

HI 

d 

-- 

H2 

= 

-r*H(z) 

0 

-U(t) 

M2 

dt 

M3 

0 

U(t) 

0 

.    M3 
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the  spins  outside  of  the  slice.  This  requirement  can  be 
expressed  by  a  desired  magnetization  profile.  For  the 
selective  180  degree  inversion  rf  pulses,  the  desired 
magnetization  can  be  expressed  as 

;;o    w2  <  i*dw 

D(l*dw)  =   <   -MO     -wl  <  l*dw  <  wl  [7] 

MO  l*dw  <  -v2 

where  D  (  1  *  d  w  )  varies  linearly  from  -w2  to  -wl  and  from  wl 
to  v  2  ;  dw  is  the  frequency  stepsize;  1  *  d  w  =  r * H ( z ) ,  where 
1  =  1 ,  2 ,...., L  and  L,  is  the  total  frequency  steps.  In  the 
actual  calculation,  the  choice  of  L  depends  on  the 
compromise  between  the  precision  requirement  and  the 
computer  time.  If  L  is  too  small,  for  example  L  =  32  as 
O'Donnell  and  Adams  (9)  chose,  there  will  not  be  enough 
precision  in  final  calculational  results  to  reflect  the 
physical  nature  of  the  problem  to  be  solved.  This  might  be 
the  reason  that  O'Donnell  and  Adams  could  not  find  the 
correct  results  even  though  they  were  the  first  to  use  the 
conjugate  gradient  method.  If  L  is  too  large,  the  computer 
time  will  be  too  long.  Because  one  should  try  several 
different  values  for  each  parameter  to  find  the  correct 
one  for  a  special  calculation,  if  L  is  too  large,  one  will 
lose  the  patience  to  try  several  times  and  lose  the  chance 
to  find  the  correct  results.  In  the  following  calculation 
L  is  chosen  as  159  in  the  optimization  procedure  of 
designing  the  90  degree  and  180  degree  selective  rf  pulses; 
L  is  chosen  as  1000  in  the  discussion  of  the  180  deeree 
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selective  refocusing  rf  pulses,  because  in  the  situation 
of  the  refocusing  rf  pulses  one  needs  nuch  nore  details  to 
understand  the  physical  nature  of  the  problem.  The 
difference,  w  2  -  wl,  is  the  width  of  the  edge  of  the 
desired  slice  profile.  The  thickness  of  the  slice  can  be 
measured  from  the  middle  point  of  a  slice  edge  to  the 
middle  point  of  the  other  edge  of  the  slice  and  it  can  be 
expressed  as  2*wl  +  2*(w2-wl)/2  =  w2  +  wl .  An  important 
parameter  for  the  slice  quality  can  be  defined  as  the 
ratio  of  the  width  of  the  slice  edge  to  the  thickness  of 
the  slice  (WTR)  and  it  can  be  expressed  as  V.'TR  = 
(w2-wl ) / ( w2+wl ) .  The  parameter  WTR  will  be  discussed  in 
detail  in  Chapters  IV  and  VI.  Also,  the  desired 
magnetization  can  be  expressed  graphically  as  shown  in 
Fig.  1. 

The  purpose  of  designing  a  selective  rf  pulse  is  to 
find  an  optimal  selective  rf  pulse,  K 1 ( t )  ,  that  causes  the 
magnetization  at  the  end  of  the  pulse,  M3(tf,l,U),  to 
approach  the  desired  magnetization,  D(l*dw),  as  closely  as 
possible,  where  tf  is  the  end  of  the  ISO  degree  pulse.  In 
optimal  control  theory  (13),  a  mathematical  function, 
which  is  called  the  objective  function  and  is  a  measure  of 
the  difference  between  the  desired  response  and  the  actual 
response,  must  be  formed  to  quantify  the  quality  of  the 
control.  One  can  use  the  following  objective  function  to 
measure  the  difference  between  D  (  1  *  d  w )  and  M 3 ( t f , 1 , U ) : 
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tf 


J(U)  »   (    dt  S'  [K3(t,l,U)-D(l*dw)]**2 

J  tO     1 


+  S'  Rl*[M3(tf ,1,U)-D(l*dw)3**2 

1 


[8] 


where  tO  is  the  beginning  of  the  130  degree  rf  pulse;  El 
are  some  numerical  values  to  be  chosen;  the  primed  capital 
letter,  £',  means  summation;  [K3(t,l,U)-I}(l*dw)]**2  means 
squared  error  distance.  Thus,  the  purpose  of  designing  a 
130  degree  inversion  rf  pulse  is  to  find  a  rf  pulse,  U (t), 
which  will  minimize  the  J(U(t))  under  the  restriction  of 
Eq.  [6]  with  the  initial  condition 

Ml  (tO)  =  0, 

M2  (tO)  =  0, 

H3  (tO)  =  MO. 
This  goal  can  be  reached  efficiently  by  the  conjugate 
gradient  method,  which  will  be  discussed  in  the  folloving 
Chapters  III  and  V. 


CHAPTER  III 
CONJUGATE  GRADIENT  METHOD  AND  OPTIMAL  150 
DEGREE  SELECTIVE  INVERSION  RE  PULSES 

To  make  the  objective  function,  Eq.  [8],  easier  to 

handle,  one  can  define  a  new  variable,  M4(t,l,U),  by  the 

equation 

dM4(t,l,U) 

=  [K3(t,l,U)  -  D(l*dw)]**2         T9] 

dt 

with  the  initial  condition  H4(tO,l,U)  =  0.  Substituting 

Eq.  [9]  into  the  objective  function,  Eq.  [8],  it  becomes 

(   tf   d 
J(U)  =  S'  \    — M4(t,l,U)dt  +  S*Rl*[M3(tf , 1 ,U)-D(l*dw) ]**2 

1  Ug  dt  1 

=  S' {M4(tf ,1,U)  +  Rl*[M3(tf tlfU)-B(l*dw)]**2}.    [10] 

1 

Combining  Eq.  [6]  with  Eq.  [9],  the  dynamic  equations 

become 

dMl(t,l,U) 

=  l*dw*M2(t,l,U) 

dt 

dM2(t,l,U) 

=  -l*dw*Kl(t,l,U)-U(t)*M3(t,l,U) 

dt 

dM3(t,l,U) 

=  U(t)*M2(t,lfU) 

dt 

dM4(t,l,U) 

=  [M3(t,l,U)-D(l*dw)]**2,  [11] 

dt 

which  can  be  called  the  extended  Dloch  equation.  The 
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[12] 


initial  condition  for  E q .  [11]  is 

Ml(tO)  ]    [  0 

K2(t0)        0 

M3(t0)        MO 

M4(t0)  J    L  0 
To  make  the  discussion  clearer,  the  dynamic  equations 
can  be  generally  written  as 


— Mi(t,l,U)=fi(Ml(t,l,U),M2(t,l,U),M3(t,l,U),H4(t,l,U),U) 

dt 

[13] 

where  i  =  1,  2,  3  and  4  ,  and  the  objective  function  can  be 
generally  written  as 
J(U)  = 

J(  Ml  ( tf ,1,U ),..., M4(tf,l.U) H 1 ( tf, L,U ),..., H4(tf,L,U)) 

[14] 
To  minimize  the  objective  function,  J ( U ) ,  one  can 
consider  a  family  of  selective  rf  pulses,  U(t)+nS(t), 
which  are  labelled  by  the  variable  q .  The  tine  function 
S(t)  is  a  modification  to  the  selective  rf  pulse  U ( t ) .  If 
the  minimizing  selective  rf  pulse  is  U(t),  the  rf  pulse 
which  makes  the  numerical  value  of  the  objective  function, 
J(U(t)  +  qS(t)),  as  snail  as  possible  will  accompany  the 
variable  q  =  0.  One  can  apply  the  standard  method  of  the 
elementary  calculus  to  determine  the  minimum.  A  necessary, 
but  not  sufficient,  condition  for  a  ninimum  of  a 
functional  is  the  vanishing  of  the  first  order  derivative 
of  the  functional.  Thus,  to  find  the  minimum  of  the 
objective  function,  one  lias 
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dJ(  U(t)+qS(t)  ) 


=  0 


i=0 


[15] 


Equation  [15]  suggests  that  one  should  find  an  expression 

dJ(U(t)+qS(t)) 

for first. 

dq 

Using  the  general  forn  of  the  objective  function, 

Eq .  [14],  and  the  chain  rule  of  the  calculus,  one  can  hav< 

dJ(U+qS)         L  4     dJ(U+qS)     d,"i  ( tf ,  1 ,  U+qS ) 
=  s'S' 

dq      q  =  0    1  i  d:Ii(  tf  ,  1  ,U+qS )         dq  q  =  0. 


[16] 


To  find  each  expression  of  the  rip,  ht  side  of  E  q .  [16],  one 
can  integrate  Eq.  [13]  and  obtain 


Mi(t.l.U)  -  Mi(tO,l,U)  +     fi(Ml MA,U)dv.     [17 

J  to 

Replacing  U  by  U+qS,  Eq.  [17]  becomes 


Mi(t,l,U+qS)  =  Mi(tO,l,U+qS)  +\   f i(Ml f . . . ,M4, U+qS)dv. 

>  to 

[18] 

Differentiating  E q .  [IS]  with  respect  to  q,  one  obtains 


—  Mi(t,l,U+qS) 

dq 

d  ,   t    4   dfi  dMj    dfi  d(U+qS) 

=  —  Mi(tO,l,U+qS)  +  1    [  S' + ]dv 

dq  J  tO   j   dMj   dq     dU     dq 


t    4   dfi  dlij    dfi 

[  s' + 5(v)  ]dv 

tO   j   dMj   dq     dU 


[19] 


where  —  Mi ( tO , 1 , U+qS )  =  0  because  Mi ( tO , 1 , U+qS )  are  the 
dq 
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equilibrium  magnetizations  which  are  constants. 

Differentiating  both  sides  of  Eq.  [19]  with  respect  to  t, 

one  can  finally  derive  a  set  of  differential  equations 

d  dMi(t,l,U+qS)    4   dfi  dMj(t ,1 ,U+qS)    dfi 

=  s'  + S(t),    [20] 

dt       dq  j   dMj        dq  dU 

where  i  =  1,  2  ,  3  and  4 .  By  the  matrix  form,  the  Eq.  [20] 

can  be  rewritten  as 

dfi    dfi  i  r  d:;i 

—  —  — .   #  #  #   —  —  —  BMI 

dill      dh'4       dq 


d 
dt 


dMl(t,l,U+qS) 

dq 


dM4(t,l,U+qS) 
dq 


df4 
dMl 


df4 
dM4 


dH4 

dq 


"  df  1  " 

dU 

+ 

• 

• 

df4 

.   dt!     . 

(t) 


[21] 


or  in  a  compact  form,  Eq.  [21]  can  be  written  as 


d  dM(t,l,U+qS)  dM 
=  fm(t,l,U+qS)  —  +  fu(t,l,U+qS)S(t) 

dt       dq  dq 


[22] 


where  M  has  four  components:  Ml,  M2,  M3  and  M4; 
dMl 


dM 
dq 


fm 


dq 

• 

dM4 
dq 
dfi 
dMl 

• 

df4 

dMl 


[23] 


dfi 
dM4 

df4 
dM4 


[24] 
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fu 


[25] 


df  1 

dU 

df4 
dU 


The  solution  of  Eq.  [22]  is 

dM(t,l,U+qS)  dM(tO,ltU+qS)    ,   t 

=  K(t,tO) +  \    K(t,v)fu(v)S(v)dv 

dq  dq         3  tO 


[26] 


where  K ( t , 1 0 )  is  the  n  x  n  matrix  solution  of  the  matrix 

differential  equation 

dK(t.v) 

=  fm(t)K(t,v); 

dt 

K(v,v)  =  I 
where  I  is  an  unit  matrix.  This  can  be  proved  easily  by 
differentiating  botli  sides  of  Eq.  [26]  with  respect  to  t 
and  combining  Eq.  [27]  as  the  following: 

d  d_M    dK(t,tO)  dM(tO) 

dt  dq       dt       dq 

t  dK(t,v) 

fu(v)S(v)dv  +  K(t,t)fu(t)S(t) 

tO   dt 


[27] 


-  fra(t)K(tftO)- 


dM(tO) 


dq 


fm(t)K(t,v)fu(v)S(v)dv  +  fu(t)S(t) 

to 

dM(tO)    (   t 

=  fm(t)[K(t,tO)  +  1   K(t,v)fu(v)S(v)dv] 

dq       )  tO 
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+  fu(t)S(t) 

dM(t) 

=  fm(t) +  fu(t)S(t). 

dq 

The  left  side  of  above  equation  and  the  right  side  of  the 
last  line  of  above  equation  consistute  an  equation  which 
is  identical  to  the  Eq.  [22],  Thus,  the  proposal  that 
Eq .  [26]  is  the  solution  of  Eq .  [22]  is  proved. 

Since  i !  (  1 0  , 1  ,  U  +  q  S  )  is  the  equilibrium  magnetization, 
dii(tO) 
dq 
[26]  becomes 

dM(tf ,l,U+qS) 

dq 
tf 


thus 


=  0 ;  letting  t  =  tf  and  combining  Eq.  [25], 


K  ( t  f  ,  v )  f  u  (  v  ,  1 ,  IJ  +  qS  )  S  (  v )  d  v 


to 


tf 


dv  S(v) 


tO 


tf 


dv  S(v) 


to 


K 1 1 ... K 1 4 


KA1.. .K44 


4      df  j 

S'  Klj  

j       dU 


dfl 

diJ 

■ 

df4 

L  dU 


[28] 


4       df  j 

S1  K4j  

j     "   dU 

where  Kij  is  the  ith  row,  jth  column  element  of  the  matrix 

K.  In  the  component  form,  Eq.  [2  8]  can  be  written  as 


21 


dMi(tf  ,1,'J+qS)    f  tf         4       df  j(v,l,U+qS) 

=  \    dv  S(v)  S'  Kij . 

dq  3  tO         j  dU 

Combining  Eq.  [16]  with  Eq.  [29],  one  gets 

dJ(U+qS)  dJ     dMi(l) 

=  S'  S' 

dq        1   i   dMi(l)    dq 

dJ    (   tf  df j(l) 

=  s*  S' \    dv  S(v)  S'  Kij 

1   i   drii(l)  J  tO         j         dU 


[29] 


tf              dfj(l) 
dv  S(v)  S'  S1  S' 


dJ 


■ij 


to 

tf 

to 


l  j    du   i  d;;i(i) 

dv  S(v)  G(v,U+qS) 


L   A   dfj(l,U)  k      dJ(v,U+qS) 

where  G(v,U+qS)  =  S'  S*  S'  Kij. 

1   j      dU     i     dMi(l) 


[30] 


[31] 


The  function  G(v,U+ qS)  is  called  the  gradient  of  the 
objective  function.  Equation  [30]  is  the  wanted  expression 
for  dJ(U+qS)/dq.  It  will  be  used  to  prove  that  the 
conjugate  gradient  algorithm  is  able  to  minimize  the 
objective  function,  J(U). 

With  the  conjugate  gradient  algorithm,  as  with  any 
convergence  algorithm,  the  first  element  of  the  convergent 
sequence,  U(0,t),  is  simply  guessed.  The  remaining  members 
of  the  sequence  are  derived  from  the  following  Eqs.  [32]  - 
[35],  Equations  [32]  -  [35]  are  called  the  conjujate 
gradient  algorithm  (13). 

S(0,t)  =  -G(t,L'(0,t)),  [32] 

U(i+l,t)  =  U(i,t)  +  qS(i,t),  [33] 

where  q  >  0  is  chosen  to  minimize  the  objective  function 
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J ( U ( i , t )  +  qS(i,t))  and  the  minimizing  q  can  be  written 
as  q(i),  i.e.  U(i+l,t)  =  U(i,t)  +  q(i)S(i,t), 

S(i+l,t)  =  -G(t,U(i+l,t))  +  b(i)S(i,t),  [34] 


b(i)  = 


tf 

i 

tO 

tf 


dt  G(t,U(i+l,t))**2 


>  0. 


[35] 


tO 


dt  G(t,U(i,t))**2 


One  can  use  the  method  of  reduction  to  absurdity  and 
prove  that  the  algorithm  Eqs.  [32]  -  [35]  will  make  tiie 
objective  function  smaller  and  smaller.  Eventually  the 
objective  function  will  reach  a  minimum.  One  can  start 
from  any  guessed  initial  rf  pulse  U ( i , t )  .  One  can  assume 
that  no  q  >  0  exists  for  which 

J(U(i,t)  +  qS(i,t))  <  J(U(i,t)).  [36] 

Therefore  one  must  have 

J(U(i,t)  +  qS(i,t))  >  J(U(i,t)).  [37] 

The  right  side  of  E  q  .  [32]  can  be  moved  to  the  left  side 
and  divided  by  q  >  0,  one  has 

J(U(i,t)  +  qS(itt))  -  J(U(i,t)) 

q 

riJ(U(i,t)  +  qS(i,t)) 
dq 
According  to  Eq.  [30]  and  combining  Eq .  [34],  S(i,t)  = 
-G(t,U(i,t))  +  b(i-l  )S(i-l , t)  ,  one  can  have 
dJ(U(i)  +  qS(i)) 

dq  q  =  0 


>  0. 


q  =  0 


[3S] 
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tf 

dv  S(i,v)  G(v,U(i)  +  qS(i)) 

tO  q»0 


tf 


dv  S(i,v)  G(v,U(i)) 


to 

tf 

to 

tf 


dv  (-G(v,U(i))  +  b(i-l)S(i-l,v))G(v,U(i)) 


tf 
dv  G(v,U(i))**2  +  b(i-l)  \   dv  S(i-1 , v)G(v,U(i) ) . 

tO  -1  to 

[39] 

One  can  combine  Eq.  [30]  with  Eq.  [33]  and  find  the  value 

of  the  second  term  of  the  right  side  of  Eq.  [39]  as  the 

f ol lowing : 

dJ(U(i-l)+qS(i-l)) 

q-q(i-l) 


tf 

i 

tO 

tf 

to 


dq 


dv  S(i-l,v)G(v,U(i-l)  +  qS(i-l)) 


dv  S(i-l,v)G(v,U(i)). 


q=q(i-l) 


[40] 


Since  q  is  chosen  to  minimize  J(U(i-l)  +  qS(i-l))  and  the 
minimizing  q  can  be  written  as  q(i-l),  one  has 
dJ(U(i-l)  +  qS(i-l)) 

dq  q=q(i-l) 

Therefore  the  value  of  the  second  term  of  Eq.  [39]  i 
tf 


=  0. 


[41] 


dv  S(i-l,v)G(v,U(i))  =  0, 


[42] 


tO 


Substituting  Eq .  [42]  into  Eq .  [39],  one  can  have 


dJ(U(i)  +  qS(i>) 
dq 


q  =  0 


tf 


dv  G(vfU(i))**2  <  0. 


[43] 


to 
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One  may  notice  that  E  q  .  [43]  contrr.  diets  E  q .  [38]  unless 
0(v,U(i))  =  0.  Thus  if  G(v,U(i))  is  not  equal  to  zero, 
then  conjugate  gradient  algorithm  will  make  the  objective 
function  smaller  and  smaller. 

To  apply  the  conjugate  gradient  algorithm,  one  needs 
to  find  an  expression  for  G(t,U(i))  which  can  be  valued 
from  the  objective  function  and  the  Ploch  equation  by  the 
numerical  procedure.  One  can  start  from  E  q  .  [31]  and 
define  a  new  variable 


dJ 
Qj(t)  =  S'  —  Kij 
i   dMi 


[44] 


or  by  matrix  form 


01 

02 

03 

.    Q4    . 

1:11  ...  K4i 


L  K14  ...  K44 


dJ 
dMl 


dJ 

d!14  J 

In  a  compact  form,  Eq.  [44]  can  be  written  as 

Q(t)  =  K+(tf,t)Jm 

where  +  means  transpose;  01,  Q2,  03  and  Q4  are  the 

components  of  the  vector  0(t);  the  letter  m  of  Jm  means 

differentiation  with  respect  to  Ml,  M2,  M3  and  M4 .  To 

obtain  a  soluble  equation  for  Q(t),  letting  S(t)  =  0  and 

setting  y(t)  =  df!(t)/dq,  Eq.[22]  becomes 

dy(t) 

=  fm(t)y(t). 

dt 

According  to  Eq.  [26],  the  solution  of  Eq.  [46]  is 


[45] 


[46] 
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y(t)  =  K(t,tO)y(tO),  [47] 

letting  t  =  tf,  tO  =  t,  E  q .  [47]  becomes 

y(tf)  =  K(tf,t)y(t).  [43] 

Combining  E  q .  [47]  and  E  q .  [48],  one  can  get 

y(tf)  =  K(tf ,t)X(t,tO)y(tO).  [49] 

Letting  t  =  tf,  Eq.  [47]  becomes 

y(tf)  =  K(tf,tO)y(tO).  [50] 

Comparing  Eq.  [50]  with  Eq.  [49],  one  can  conclude  that 

K(tf,tO)  =  K(tf ,t)K(t,tO).  [51] 

Assuming  that  tO  and  tf  are  constants,  differentiating 

Eq.  [51]  with  respect  to  t  and  cor.  bin  in  g  Eq.  [27],  one 

obtains 

dK(tf,t)  dK(t.tO) 

0  = K(t,tO)  +  K(tf,t) 

dt  dt 

dK(tf.t) 

= K(t.tO)  +  K(tf,t)fr.(t)K(t,tO) 

dt 

d?:(tf,t) 

=  ( +  K(tf,t)fm(t))K(t,tO).  [52] 

dt 

Since  K ( t , 1 0 )  is  not  equal  to  zero,  one  gets  from  Eq.  [52] 

dX(tf ,t) 

=  _K(tf ,t)fm(t).  [53] 

dt 

According  to  Eq.  [27],  the  terminal  condition  of  Eq.  [53] 

is 

K(tf,tf)  =  I.  [54] 

Taking  the  transpose  of  Eq.  [53]  and  Eq.  [54]  produces 

dK+(tf,t) 

=  -fm+K+(tf,t); 

dt 
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K+(tf,tf)  =  I. 


[55] 


Postmultiplying  both  sides  of  En.  [55]  by  Jm  gives 

dK+(tf,t) 

jn,  =  -fm+K+Ctf  ,t)Jm; 

dt 

K+(tf,tf)Jm  =  Jn.  [56] 

By  the  definition  of  Q(t),  Eq.  [45],  Eq.  [56]  becomes 

dQ(t) 

=  -fra+  Q( t ) ; 

dt 

0(tf)  =  Jm.  [57] 

Equation  [57]  is  the  soluble  equation  for  C(t)  and  is 

called  an  adjoint  equation  of  the  Eq.  [27]. 

After  solving  Eq.  [57],  one  can  finally  get  the 

gradient  by  Eq.  [31] 


dfj(l.U) 
G(t,U)  =  S'  s'  Qj(t). 

1   j     dU 


[53] 


-fn  +  = 


Considering  the  extended  Eloch  equation,  Eq.  [11], 
the  definition  equations,  Eq.  [24]  and  Eq.  [25],  become 
0     l*dw      0  0 

-l*dw    0     -U(t)  0 

0     U(t)     0     -2(M3(t)-D(l*dw)) 
0      0       0  0 

0 
-H3(t) 
M2(t) 
0 

Considering  the  objective  function,  Eq.  [10],  the 
terminal  condition  of  Eq.  [57]  becomes 


[59] 


fu(t)  = 


[60] 
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Q(tf)  = 


0 


0 


:El*(M2(tf)-P(l*dw)) 


] 


[61] 


According  to  E  q  .  [60]  and  Eq.  [61],  the  gradient  of 

the  objective  function,  Eq.  [Si1],  becones 

L 
G(t,U(i))  =  S'  (-M3(t,l)Q2(t,l)  +  M2(t,l)Q3(t,l)).     [62] 


In  summary,  the  procedure  for  designing  the  ISO 
degree  selective  inversion  rf  pulses  is  in  the  following : 

1.  Guessing  an  initial  180  degree  rf  pulse,  for  ex an pie, 
a  150  degree  truncated  sine  function  rf  pulse,  integrate 
trie  Eq.  [11]  under  the  initial  condition  Eq.  [12]  and  find 
the  first  value  of  the  objective  function,  J(U,0). 

2.  According  to  the  solution  of  Eq.  [11]  and  using  the 
Eq.  [59]  and  Eq.  [61],  construct  the  adjoint  equation 

Eq.  [57]  and  integrate  Eq.  [57]  backwards  in  time  from  tf 
to  tO  to  get  the  vector  function  0(t). 

3.  Construct  the  gradient  of  the  objective  function  by 
Eq.  [621. 

4.  Using  Eq.  [32]  and  Eq.  [23],  construct  a  new  rf  pulse, 
U(t).  The  numerical  coefficient,  q,  is  chosen  to  minimize 
the  objective  function  J  for  the  given  U(i,t)  and  S(i,t) 
and  the  minimizing  q  can  be  found  by  using  the  Golden 
Section  Search  (16).  The  new  rf  pulse  is  use:!  to  find  a 
new  value  of  the  objective  function,  J(n).  Comparing  the 
new  value  with  the  old  value  of  the  objective  function, 
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J(n-l),  stop  here  if  J(n)  >  J(n-l),  or  construct  a  new 
2  r  a  c!  i  n  e  t  and  continue  by  u  s  i  n  g  E  q .  [34]  if  J  (  n  )  <  J  (  n  - 1 )  , 
To  solve  the  extended  i:loch  equation,  Eq.  [11],  and 
the  adjoint  equation,  Eq.  [52],  numerically,  the  fourth 
or (sr  Runge-Kutta  method  can  be  used  which  is  in  the 
Appendix  A .  The  four  steps  above  can  be  implemented  by  a 
Fortran  program  that  is  in  the  Appendix  B.  The  parameter 
values  used  in  the  program  of  Appendix  B  to  obtain  the 
optimal  rf  pulses  in  Fig.  2,  Fig.  11  and  Fig.  14  of 
Chapter  IV  are  in  a),  b)  and  c)  of  Appendix  C. 


CHAPTER  IV 
EXPERIMENTAL  STUDY  OF  THE  OPTIMAL  180 
DEGREE  SELECTIVE  RF  PULSES 

A  detailed  derivation  and  procedure  of  how  to  use  the 
conjugate  gradient  method  to  design  the  180  degree 
selective  inversion  rf  pulses  is  already  given  in  the 
Chapter  III.  In  this  chapter  several  special  cases  of  180 
degree  selective  rf  pulse  will  be  discussed  and  the 
simulated  results  will  be  compared  with  the  experimental 
results.  The  experimental  results  have  confirmed  the 
simulation  prediction  that  the  optimal  180  degree 
selective  rf  pulses  have  much  better  selectivity  than  the 
widely  accepted  sine  function  rf  pulses  in  the  both 
situations  of  longitudinal  magnetization  inversion  and 
transverse  magnetization  refocusing.  An  important 
conclusion  of  this  chapter  is  that  the  responses  to  a 
selective  180  degree  rf  pulse  will  be  the  same  no  matter 
how  it  is  used:  as  a  selective  inversion  pulse  or  as  a 
selective  refocusing  pulse.  Thus  one  can  apply  the  design 
procedure  for  the  180  degree  selective  inversion  rf  pulse 
to  obtain  the  selective  rf  pulse  for  the  transverse 
magnetization  refocusing. 

Before  starting  a  computation,  one  should  choose  a 
desired  slice  profile  and  an  initial  trial  function.  A 
desired  slice  profile  for  a  180  degree  selective  inversion 
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rf  pulse  is  already  shown  in  Fig.  1,  The  slice  thickness 
is  measured  in  frequency  unit.  As  it  was  already  pointed 
out  in  the  Chapter  II,  the  important  parameter  of  the 
slice  quality  is  the  ratio  of  the  width  of  the  slice  edge 
to  the  thickness  of  the  slice  (WTR)  and  it  can  be 
represented  by  WTR  =  ( w2-wl ) / ( wl+w2 ) ,  where  w2-wl  is  the 
width  of  the  slice  edge  and  wl+w2,  which  is  equal  to  2*v;l 
+  2*(w2-wl)/2,  is  the  half  height  thickness  of  the  slice 
measured  from  the  middle  point  of  an  edge  of  the  slice  to 
the  middle  point  of  the  other  edge.  It  is  obvious  that  the 
smaller  the  WTR  value,  the  better  the  slice  quality,  if 
all  other  conditions  are  the  same.  A  four  zero  truncated 
sine  function,  as  shown  in  Fig.  2  (dotted  line),  is  chosen 
as  an  initial  trial  function  for  the  conjugate  gradient 
algorithm.  The  inversion  response  to  the  180  degree  four 
zero  sine  function  rf  pulse  is  shown  in  Fig.  3  (dotted 
line).  If  the  slice  thickness  of  a  desired  slice  profile 
is  identical  with  that  of  the  response  to  the  four  zero 
sine  function  pulse  and  if  one  chooses  WTR  =  2/9,  as  shown 
in  Fig.  1,  the  corresponding  optimal  inversion  rf  pulse  is 
shown  in  Fig.  2  (solid  line)  and  the  inversion  response  to 
the  optimal  rf  pulse  is  shown  in  Fig.  3  (solid  line). 

Since  a  selective  inversion  rf  pulse  will  be  used  as  a 
selective  refocusing  rf  pulse  in  the  following  discussion, 
the  solution  of  the  Bloch  equation  without  relaxation  in 
the  refocusing  situation  will  be  considered  here.  For 
explanation,  one  can  consider  a  simple  situation  in  which 
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there  is  no  selective  rf  pulse  and  there  is  only  a 
z-direction  magnetic  field  gradient.  In  solving  the  Bloch 
equation  without  relaxation  in  the  rotating  frame  by 
computer,  there  are  1000  points  in  the  gradient  dimension 
and  each  point  corresponds  to  a  solution  of  the  Bloch 
equation  without  relaxation  in  which  the  specific 
z-direction  magnetic  field  is  determined  by  the  gradient. 
If  the  same  initial  transverse  magnetization  along  the 
y-direction  is  used  for  each  solution,  the  y-components  of 
the  solutions  for  1000  points  can  be  plotted  together  and 
are  shown  in  Fig.  4.  The  dephasing  action  of  the  gradient 
is  not  included  in  this  simulation.  If  the  optimal 
selective  rf  pulse  of  Fig.  2  is  used  as  a  selective 
refocusing  pulse  for  the  transverse  magnetization,  the 
response  to  the  optimal  pulse  is  shown  in  Fig.  5.  For 
comparison,  the  solid  line  of  Fig.  3  and  the  line  of 
Fig.  5  are  plotted  together  in  Fig.  6.  Two  lines  of  Fig.  6 
fit  together  very  well.  If  the  sine  function  rf  pulse  of 
Fig.  2  is  used  as  the  selective  refocusing  pulse  for  the 
transverse  magnetization,  the  response  to  the  sine 
function  pulse  is  shown  in  Fig.  7  (dotted  line).  The  solid 
line  of  Fig.  7  is  the  dotted  line  of  Fig.  3.  Two  lines  of 
Fig.  7  fit  together  very  well  also.  Thus,  the  responses  to 
a  selective  180  degree  rf  pulse  will  be  the  same  in  either 
situation  of  use:  as  a  refocusing  pulse  or  as  an  inversion 
pulse.  In  the  following,  the  selective  ISO  degree  rf 
pulses  will  be  used  experimentally  as  transverse 
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magnetization  refocusing  pulses  and  the  experimental 
results  will  be  compared  with  the  simulation  results  when 
the  180  degree  pulses  are  used  as  longitudinal 
magnetization  inversion  pulses.  The  Fortran  program  for 
selective  refocusing  rf  pulse  study  is  in  Appendix  D. 
The  instrument  used  in  the  experiments  is  the  CE 
CSI-II  2T  spectrometer /imager  system.  The  sample  is  a  tube 
of  water  and  it  is  along  the  main  magnetic  field 
( z-direction) .  The  tube  length  is  100  mm  and  the  tube 
diameter  is  10  mm.  A  spin  echo  pulse  sequence,  as  shown  in 
Fig.  8 ,  is  used  in  the  experiments.  In  the  pulse  sequence, 
the  90  degree  rf  pulse  is  not  selective  and  180  degree 
refocusing  rf  pulse  is  frequency  selective.  The  pulse 
sequence  in  Fig.  8  results  in  a  one-dimensional  projection 
image  in  the  z-direction.  In  the  experiments,  the 
non-selective  90  degree  pulse  length  is  165  us,  the  four 
zero  sine  function  selective  180  degree  pulse  length  is 
4.5  ms  and  the  six  zero  sine  function  selective  180  degree 
pulse  length  is  4.5  ms  as  well.  The  selective  180  degree 
rf  pulse  lengths  can  be  reduced  by  increasing  the  output 
power  level  of  the  rf  pulse  amplifier.  The  echo  time,  TE, 
is  30  ms  and  the  acquisition  data  block  size  is  1024 
points.  The  experimental  results  for  the  180  degree  four 
zero  sine  function  rf  pulse  and  the  corresponding  optimal 
rf  pulse  are  shown  in  Fig.  9  and  Fig.  10,  respectively. 
Comparing  Fig.  9  and  Fig.  10  with  Fig.  3,  one  can  see  that 
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the  experimental  results  are  very  close  to  the  computer 
simulations . 

The  reason  why  there  is  not  any  signal  in  the 
experimental  results  of  Fig.  9  and  Fig.  10  that 
corresponds  to  the  oscillating  magnetization  of  Fig.  6  and 
Fig.  7  is  the  gradient  dephasing  action.  Applying  a  pulse 
sequence  in  which  there  is  no  180  degree  selective 
refocusing  pulse  and  the  rest  is  the  same  as  in  Fig.  8, 
one  cannot  get  any  signal  experimentally,  because  the 
applied  gradient  of  the  pulse  sequence  will  dephase  the 
magnetization.  By  solving  the  Rloch  equation  without 
relaxation,  the  y-component  of  the  magnetization  for  the 
pulse  sequence  without  selective  rf  pulse  is  already  shown 
in  Fig.  4 .  Thus,  Fig.  4  corresponds  to  zero  signal 
obtained  experimentally  by  the  pulse. sequence  without 
selective  pulse.  By  applying  the  pulse  sequence  of  Fig.  8, 
the  selective  refocusing  180  degree  pulse  will  select  only 
a  portion  of  the  oscillating  magnetization  of  Fig.  4.  The 
refocusing  gradient  has  the  refocusing  action  only  to  this 
portion  of  the  magnetization.  The  magnetization  that  is 
not  selected  by  the  180  degree  pulse  will  be  dephased. 
Thus,  one  can  only  obtain  the  signal  from  the  selected  and 
refocussed  magnetization. 

If  the  thickness  of  a  desired  slice  is  identical  with 
that  of  the  response  to  the  four  zero  sine  function  rf 
pulse  and  if  one  reduces  WTR,  that  is,  one  tries  to 
produce  a  sharper  edge,  for  example,  one  chooses  WTR  =  1/7, 
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the  optimal  rf  pulse  obtained  is  shown  in  Fig.  11  (dotted 
line).  The  corresponding  response  to  the  dotted  line  pulse 
of  Fig.  11  is  shown  in  Fig.  12  (dotted  line).  One  can 
notice  that  even  though  there  are  sharper  edges  on  the 
slice,  there  are  more  artifacts  on  the  top  of  the  slice 
and  outside  of  the  dice.  Thus  one  may  conclude  that  for  a 
particular  pulse  length  there  is  a  minimum  value  of  WTR, 
represented  by  WTRmin,  at  a  certain  slice  thickness.  If 
WTR  >  WTRmin,  one  can  always  find  an  rf  pulse  to  which  the 
response  will  be  very  close  to  the  desired  slice.  If  WTR  < 
WTRmin,  one  can  obtain  a  sharper  slice  edge,  but  the  price 
is  that  there  are  more  artifacts  on  the  top  of  the  slice 
and  outside  of  the  slice,  as  shown  in  Fig.  12.  The  above 
observation  has  been  proved  by  experimental  results  as 
shown  in  Fig.  10  and  Fig.  13.  If  one  changes  the  pulse 
length  while  retaining  the  shape,  the  slice  thickness  will 
be  changed.  For  example,  if  the  pulse  length  of  a  IPO 
degree  four  zero  sine  function  rf  pulse  is  doubled  and  the 
pulse  shape  is  still  four  zero  sine  function,  the  slice 
thickness  will  be  half,  if  all  other  conditions  are  the 
same.  However,  the  value  of  WTR,  which  determines  the 
sharpness  of  a  slice,  does  not  depend  on  the  pulse  length 
and  will  not  be  changed  if  only  the  pulse  length  is 
changed.  From  the  discussion  above  one  may  conclude  that 
one  cannot  make  a  slice  as  sharp  as  one  wants  by  a  single 
selective  rf  pulse  at. a  certain  slice  thickness.  The 
property  from  the  discussion  above  can  be  attributed  to 


44 


in 

* 

4-1 

If4 

(U 

c 

t- 

u 

CD 

•H 

u 

i— 1 

r-J 

0J 

cd 

cc 

ai 

CO 

X! 

4-1 

E 

c 

4-> 

VM 

•H 

•H 

•o 

■H 

4J 

t-H 

oj 

c 

T) 

Q. 

u 

4J 

c 

"O 

•H 

x: 

a> 

00 

CD 

4-1 

B 

4J 

CU 

CD 

•H 

ra 

4J 

-c 

C 

s 

c 

o 

CO 

c 

0J 

c. 

CO 

4J 

jr 

CO 

0) 

C 

• 

4J 

0) 

en 

v 

cs 

kl 

■-I 

DO 

c 

3 

1 

• 

a; 

SJ 

D. 

u 

OC  X! 

XI 

c 

•H 

s 

4J 

y-i 

cj 

Pn 

ki 

u 

a. 

<4-l 

<4-l 

do 

o 

.— i 

CD 

c 

1— 1 

to 

E 

TO 

c 

3 

4-1 

E 

•H 

4-1 

a. 

co 

•H 

i— I 

a 

x: 

^^ 

4-1 

X! 

i«-i 

4J 

• 
E 

Q.-U 

4-) 

ki 

o 

•H 

x: 

i— 1 

00 

^H 

4->         • 

H 

OJ 

O 

IS 

TO 

•h  r- 

(1) 

C/j 

E 

S  \ 

k. 

0) 

-H 

1— 1 

ec 

• 

E 

4J 

1— 1 

01 

OJ 

cs 

CL 

ca    II 

■a 

i— 1 
■H 

n 

c 

u 

•H   « 

O    "4-1 

OJ 

ki 

4J    H 

00 

o 

jr 

0) 

C  3 

I— 1 

kl 

4-1 

XI 

0) 

Q. 

4.) 

-a  4-i 

o 

CO 

o 

•H     3 

> 

SJ 

•H 

c 

XI 

H 

u 

CC 

CO 

•H 

XI 

•H       « 

• 

iH 

u 

n 

0J 

1—1 

CO 

•H 

4J 

co   n 

i— 1 

-C 

C 

CO    -H 

-G 

3 

o 

CD   i-H 

• 

CD 

CO 

c 

cc 

kl 

u 

OJ 

.*  -3 

•H 

•H 

CO 

k< 

O    -H 

Pn 

09 

T— 1 

arid 

a; 

3 

CU 

X!    O 

•a 

a. 

u 

4J      CC 

epnj!|duuv     eAijeisy 


45 


—    o 
o 

o 
o 


o 
o 
o 
o 


0) 

U-t 

c 

o 

•H 

r-l 

oo 

V 

-a 

oo 

•H 

!-H 

iH 

a 

O 

D. 

oo 

tM 

ai 

>H 

x; 

4-1 

f-H 

CO 

o 

E 

4-> 

^^ 

•H 

N 

4J 

00 

X 

CUTS 

w 

O 

c 
o 

>. 

u 

o. 

o 

c 

9 

ja 

00 

j_) 

CD 

3 

Ih 

CT 

0 

u 

9 

4-1 

o 

*- 

(J 

U. 

00 

0) 

0) 

00 

c 

c 

•H 

o 

1—1 

D, 

00 

•o 

0) 

•H 

u 

i— 1 

o 

CD 

00 

Xi 

H 

V 

• 

H 

CN 

i-H 

•     • 
i— i  >— i 

• 

I— 1    .—1 

oo 

•H 

•        • 

&-, 

QO    O0 

2 

l 


£,  &-. 


46 


o 
o 
o 
o 


-    o 


o 

o 
o 
o 


d     ai 

I-    -C 

ai  cc  <4-i 

co  -i 

x:  -h  o 

• 

4-1 

cu  e  o 

4->    4J 

U    O  i-i 

O  H    1C 

CD 

4->   3  x: 

x:   co     • 

CO   4-> 

4J    4J     00 

CU     01 

CJ    -H 

co  u  co 

x:   co  tt.. 

e        cc 

CC4-I 

O    CU 

3    -H    <4-l 

D.JC    C 

O    4J     O 

WHO 

XI    u 

CU           -H 

4J     CO    4-1 

U        •    4-1 

CC 

r-  -h 

c   cu  x: 

ai  ■^-•c 

0)      (-1     4-1 

JS  —>   c 

>  o 

4-1             O 

tJ    E    C 

II      CJ 

co 

<4-l 

•  cu  x: 

N 

X 

O   KH 

CJ      U     4J 

E-h    CO 

Du  CO 

*^ 

u2  iJ 

03           0) 

i-H             C 

X    ID    U 

>, 

3  r  (u 

CO    1-   -H 

CJ 

K     4-1      E 

CU  i—l 

c 

CU   i-H  -H 

cu  xz   co 

0 

U     3     U 

CO     4J 

3 

CU 

iH           CU 

a 

i— 1  •-<    c 

3     « x: 

• 

CO    -H      X 

c  d'  4-i 

4-1           CD 

u 

C      • 

<W    'H«H 

CU    00  CU 

t*  iH    O 

E  -H    E 

CO 

■H  Pt-i    CO 

CU           CU 

U           CO 

x:  cu  -a 

ai  in 

4J  x:  i-i 

a.  o  ai 

4J     CO 

X          XI 

U             4-1 

0)     01    4_> 

O  «*-i    3 

CO 

"WOO 

QJH    U 

x:    3   cu 

4J   w  -a 

H     O.T3 

am  c 

c 

CU    00  CO 

•  U->    3 

CJ  -a 

cn  u 

x   cu   c 

•— i       -a 

cu        o 

t-H    0) 

1-4     4-1 

•    CO    c 

o>    CU 

Ot    E-H 

a.  cu 

•rl    -H     CC 

•     l-c      O 

t-     4-1     4J 

CC    CO    -H 

Q.X) 

•H   X:  r-l 

O    O 

En    CO    CO 

47 


the  nature  of  the  spin  system.  By  the  terminology  of 
control  theory,  one  properly  could  conclude  that  not  every 
desired  slice  profile  is  reachable  by  a  practical  rf  pulse. 

If  one  tries  to  improve  the  slice  quality,  that  is, 
to  reduce  the  value  of  WTR,  it  can  be  noticed  that  one 
must  choose  a  desired  slice  that  has  a  broader  slice 
thickness  if  the  pulse  length  is  fixed.  Because  the 
conjugate  gradient  method,  like  most  other  iterative 
numerical  procedure,  can  only  find  a  local  minimum,  one 
can  start  the  computation  from  a  trial  pulse  shape  to 
which  the  response  will  have  a  nearly  equal  slice 
thickness  to  that  of  the  desired  slice.  One  can  choose  a 
sine  function  shaped  pulse  with  more  zeros  as  the  initial 
trial  function  because  the  pulse  of  sine  function  with 
more  zeros  will  have  a  broader  frequency  response  for  a 
fixed  pulse  length.  For  example,  one  can  choose  a  desired 
slice  that  has  the  slice  thickness  which  is  identical  with 
that  of  the  response  to  a  six  zero  sine  function  rf  pulse 
and  WTR  =  1/7,  the  corresponding  optimal  rf  pulse  is  shown 
in  Fig.  14  (solid  line).  The  simulated  response  to  the 
optimal  pulse  is  shown  in  Fig.  15.  The  experimental 
results  are  shown  in  Fig.  16  and  Fig.  17,  respectively. 
The  slice  edge  in  Fig.  17  is  sharper  than  that  in  Fig.  10; 
however,  the  slice  in  Fig.  17  has  a  broader  frequency 
range.  If  one  wants  the  same  spatial  slice  thicknesses  for 
two  slices  at  a  specific  value  of  WTR,  the  broader 
frequency  range  requires  a  higher  gradient.  One  can  notice 
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that  the  optimal  rf  pulse  corresponding  to  the  180  degree 
six  zero  sine  function  pulse  has  a  higher  power  level  than 
the  optimal  rf  pulse  corresponding  to  the  180  degree  four 
zero  sine  function  pulse  if  the  pulse  lengths  are 
identical,  or  the  optimal  rf  pulse  corresponding  to  the 
180  degree  six  zero  sine  function  pulse  has  a  longer  pulse 
length  than  the  optimal  rf  pulse  corresponding  to  the  180 
degree  four  zero  sine  function  pulse  if  the  power  levels 
are  identical.  Thus,  one  may  conclude  that  a  higher  slice 
quality  requires  a  higher  rf  pulse  power  level  in 
addition  to  a  higher  gradient  or  a  longer  rf  pulse  length. 

Even  though  the  experimental  results  in  Fig.  10  and 
Fig.  17  are  very  close  to  the  calculation  results,  there 
are  still  some  differences  between  the  experiments  and 
calculations.  The  reason  for  the  difference  between  the 
simulations  and  the  experiments  is  in  the  following:  the 
rf  amplifier  in  our  instrument  is  not  quite  linear, 
especially  in  the  regions  near  zero  and  maximum  power 
level.  Even  though  one  can  make  some  compensations  for  the 
nonlinear i ties ,  it  is  very  hard  to  make  the  rf  pulses  in 
the  experiments  exactly  as  shown  in  Figs.  2,  11  and  14. 
Because  of  the  nonlinearity  of  the  rf  power  amplifier,  the 
optimal  rf  pulse  corresponding  to  a  more  zero  sine 
function  pulse  will  be  more  difficult  to  ;aake  since  it  has 
a  higher  power  level  and  a  more  complicated  structure  near 
the  zero  power.  The  second  reason  for  the  difference  is 
that  in  our  experiments  one  must  input  the  rf  pulse 
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amplitude  point  by  point;  the  pulse  sequence  compiler  will 
allow  at  most  about  90  points  to  make  an  experimental  rf 
pulse.  If  there  are  more  points,  the  CSI-II  system  will 
not  compile  the  pulse  sequence.  However,  90  points  are  not 
enough  to  accurately  make  a  experimental  rf  pulse.  If 
there  are  some  improvements  in  above  two  situations,  the 
experimental  results  will  be  closer  to  the  calculations. 

The  author  was  required  to  input  the  experimental  rf 
pulses  point  by  point;  however,  a  Fourier  analysis  of  the 
calculated  optimal  rf  pulses  might  be  useful  for  some 
readers  in  defining  the  optimal  rf  pulses.  The  Fourier 
coefficients  for  the  Fourier  series  representation  of  the 
optimal  rf  pulses  in  Fig.  2  and  Fig.  14  are  shown  in  the 
Appendices  E  and  F,  respectively.  The  Fortran  program  to 
obtain  the  Fourier  coefficients  of  the  Fourier  series 
representation  of  the  optimal  rf  pulses  is  in  the 
Append  ix  G . 

The  relation  between  the  slice  thickness  and  the 
gradient  magnitude  will  be  discussed  in  this  paragraph.  If 
the  pulse  sequence  of  Fig.  8  is  applied,  all  three 
intervals  of  the  z-direction  magnetic  field  gradient  will 
relate  to  the  thickness  of  the  slice.  Since  there  is  a 
quantitative  relation  between  the  dephasing  interval  and 
the  refocusing  interval,  there  are  only  two  independent 
gradient  magnitudes  that  will  determine  the  slice 
thickness.  The  sample  is  a  tube  of  water  with  a  partition 
at  the  middle  of  the  tube.  The  length  of  the  tube  is  100 
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ram  and  the  diameter  of  the  tube  is  10  mm.  The  thickness  of 
the  partition  is  about  1  mm.  The  experimental  results  is 
shown  in  Fig.  18.  If  one  doubles  the  selective  gradient 
magnitude  while  keeping  other  conditions  unchanged,  the 
slice  thickness  will  be  half  with  the  same  value  of  WTR , 
or  vice  versa,  as  shown  in  a),  b)  and  c),  d)  of  Fig.  18. 
From  a)  to  b)  of  Fig.  IS,  the  selective  gradient  magnitude 
is  doubled,  the  frequency  thickness  and  the  spatial 
thicicness  of  the  slice  are  half  simutaneotisly.  If  one 
halves  the  refoc using  gradient  magnitude  while  keeping 
other  conditions  unchanged,  the  frequency  slice  thickness 
will  be  half  with  the  same  value  of  WTR,  or  vice  versa,  as 
shown  in  a)  and  c)  of  Fig.  18.  In  this  case  it  can  be 
emphasized  that  the  spatial  thickness  of  the  slice  is  not 
changed  even  though  the  frequency  thickness  of  the  slice 
is  half.  If  one  doubles  the  refocusing  gradient  magnitude 
and  the  selective  gradient  magnitude  simultaneously,  the 
frequency  slice  thickness  will  not  be  changed  and  the 
spatial  slice  thickness  will  be  half.  Thus,  if  one  doubles 
the  selective  gradient  magnitude,  the  spatial  slice 
thickness  is  always  half,  no  natter  the  refocusing 
gradient  magnitude  is  changed  or  not.  It  can  be  clear  from 
the  Fig.  18  that  the  following  relation: 

dz  =  dv/(r*Gz)  [63] 

is  true  only  when  Gz  is  the  refocusing  gradient  and  it  is 
not  true  when  Gz  is  the  selective  gradient,  where  dz  is 
the  spatial  thickness  of  the  slice. 
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The  the  influence  of  the  gradient  nag nitudes  on  the 
signal  magnitude  will  he  discussed  in  this  paragraph.  If 
one  increases  the  selective  gradient  magnitude,  the  signal 
magnitude  will  not  be  changed  as  shown  in  a)  and  b)  of 
Fig.  IT.  If  one  halves  the  refocusing  gradient  magnitude, 
the  signal  magnitude  will  be  double  as  shown  in  a)  and  c) 
of  Fig.  18.  Thus,  only  the  refocusing  gradient  magnitude 
'..'ill  have  an  effect  on  the  signal  magnitude. 


CHAPTER  V 
CONJUGATE  GRADIENT  METHOD  AND  THE  OPTIMAL 
90  DEGREE  SELECTIVE  RF  PULSES 

As  we  will  see  later,  a  truncated  90  degree  sine 

function  rf  pulse  is  already  a  good  90  degree  selective  rf 

pulse,  but  one  can  still  use  the  optimal  control  method 

and  improve  the  rf  pulse  shape  further  to  get  better  90 

degree  slice  selectivity.  Since  the  90  degree  selective  rf 

pulse  design  incorporates  a  refocusing  gradient  after  the 

90  degree  rf  pulse,  the  problem  to  design  the  90  decree 

selective  rf  pulse  is  a  little  more  complicated  than  that 

of  the  18  0  degree  pulse  design.  The  time  interval  when 

there  is  a  90  degree  selective  rf  pulse  can  be  called  the 

selective  period.  The  interval  when  there  is  no  90  decree 

selective  rf  pulse  and  there  is  only  refocusing  magnetic 

field  gradient  can  be  called  the  refocusing  period.  As  in 

the  situation  of  the  180  degree  selective  rf  pulse  design 

one  should  decide  what  is  the  desired  90  degree 

magnetization,  and  then  derive  a  procedure  which  can  be 

programmed  to  find  the  optimal  selective  90  degree  rf 

pulse  that  causes  the  transverse  magnetization  at  the  end 

of  the  refocusing  period  to  approach  the  desired 

magnetization  as  closely  as  possible.  In  the  following 

discussion,  we  set  U(t)  =  -r*Hl(t)  as  in  the  case  of  the 

design  of  the  180  degree  selective  rf  pulse,  in  the 
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refocusing  period  U  ( t )  =  0.  Ml,  f  1 2  and  1 !  3  represent  x,  y 
and  z  components  of  the  magnetization  in  the  selective 
period,  and  Ml,  M2  and  ii'3  represent  x,  y  and  z  components 
of  the  magnetization  in  refocusing  period,  respectively. 
The  bar,  -,  does  not  mean  a  vector  here. 

Because  there  is  no  rf  pulse  in  the  refocusing  period, 
the  objective  function,  which  is  the  measure  of  the 
difference  between  the  desired  transverse  magnetization 
and  the  actual  transverse  magnetization,  can  be  expressed 
as 

f  tf 
J  =  \     dt  S'  [Ml(t,l,0)**2  +  !'.2(t,l,0)-:*2-P(l*dv)**2]**2 

)   tl     1 

+  S'  Rl*[Ml(tf ,1,0)**2  +  M2(tf ,1,0)**2-D(l*dw)**2]**2 

1 

T64] 
where  tl  and  tf  are  the  beginning  and  the  end  of  the 
refocusing  period,  D ( 1 * d w )  is  the  desired  transverse 
magnetization  and  can  be  expressed  as 

0     w2  <  l*dw 

MO   -wl  <  l*dw  <  wl  [65] 

I   0  l*dw  <  -w2 

where  w2-wl  is  the  width  of  the  slice  edge  and  wl+w2  is 
the  thickness  of  the  slice.  The  desired  transverse 
magnetization,  Eq .  [65],  can  be  graphically  shown  in 
Fig.  19.  The  length  of  the  refocusing  magnetic  field 
gradient  can  be  chosen  to  be  identical  with  the  90  degree 
selective  rf  pulse  length  for  convenience  in  the  analysis. 


L(l*dw)  = 
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Thus,  the  magnitude  of  the  refocusing  magnetic  field 
gradient  is  half  of  the  selective  magnetic  field  gradient. 
Considering  the  problem  in  the  rotating  frame  and  assuming 
the  selective  rf  pulse  is  applied  along  the  x  direction  of 
the  rotating  frame  on  resonance,  the  I' loch  equation  for 
the  selective  period  is 


d 
dt 


Ml 


l  m: 


0 


-I! 


Ml 
M2 
M3 


0      l*dw 
-l*dw     0 

0      u 

with  the  initial  condition 

M1(0)  =  0, 

K2(0)  =  0, 

M3(0)  =  MO. 
The  Eloch  equation  for  the  refocusing  period  is 


[66] 


[67] 


d 
dt 


Ml 


Ml 
M2 


[6R] 


[69] 


0       -0.5*l*dw    0 
0.5*l*dw      0       0 

0  0       0  J 

The  initial  condition  for  Eq.  [6S]  is 

Ml(tl)  =  Ml(tl), 

M2(tl)  =  M2(tl), 

M3(tl)  =  M3(tl). 
One  can  define  a  new  variable,  M  4  (  t ,  1  ,  0  )  ,  by  the  equation 

dM4(t,l,0) 

=  [Ml(t,l,0)**2  +  M2(t,l,0)**2  -  D(l*dw)**2]**2 

dt 

[70] 
with  the  initial  condition  M4(t 1,1,0)  =  0.  Substituting 
Eq.  [70]  into  Eq.  [64],  the  objective  function  becomes 
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J  = 

S1  (M4(tf,l,0)  +  R1*[M1**2  +  K2**2  -  D(l*dw)**2]**2] 

1 

Combining  Eq.  [  6 c]  ]  with  E  q .  [7  0],  the  extended  Bloch 
eauations  for  the  refocusinp  oeriod  is 


dMl 

„ 

-0.5*l*dw*M2 

dt 

dM2 

0.5*l*dw*Kl 

dt 

d7l3 



= 

0 

dt 

dM4 

=  [  M 1  *  #  2  +  M2  **  2  -  D(l*dw)**2]**2 

dt 

with  the  initial  condition 

Ml(tl)  =  Hl(tl), 

~2(tl)  =  M2(tl), 

H3(tl)  =  M3(tl), 

M4(tl)  =  0. 
Two  cascaded  systems,  Eq.  [66]  and  Eq.  [72]  can  be 
generalized  to  Eq.  [74]  and  Eq .  [75],  respectively. 
dMi(t,l,U) 

dt 
where  i  =  1,  2  and  3,  and 


dMj(t,l,0) 

=  fj(Ml(t,l,0)§...,M4(tt1.0)) 

dt 


[71] 


[72] 


[73] 


=  fi(Ml(t,l,U),;i2(t,l,U),M3(t,l,U),U)      [74] 


[75] 


where  j  =  1,  2,  3  and. 4.  Because  the  variables  in  Eq.  [75] 
are  related  to  the  variables  in  Lq.  [74]  by  the  initial 
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condition,  E q .  [73],  the  objective  function,  E  q  .  [711,  can 
be  generalized  to 

J(U)  =  J(Ml[tf,l,Ml(tl,l,U),...,H3(tl,l,U)] 

M4[tf  ,l,i:i(tl,l,U) K3(tl,l,U)])  [76] 

where  1  =  1  ,  2  ,  .  .  .  ,  L .  L  is  chosen  to  be  equal  to  15?. 

One  can  start  to  derive  the  optimization  procedure 
from  the  partial  derivative  of  the  objective  function, 
Eq.  [76],  as  we  did  in  deriving  the  optimization  procedure 
for  designing  the  180  degree  selective  rf  pulse  from  the 
E  q .  [14],  One  obtains  from  E  q  .  [76] 

dJ(U+qS)    L   4     dj      3   dMj(tf)  dMi( tl ,1 ,U+qS) 

=  s'  S1 S' . 

dq       1   j   dMj(tf)   i   dMi(tl)        dq 

[77] 

Integrating  Eq.  [75],  one  gets 

Mj(t)  -  Mj(tl)  +  (    fj(~l(v,l),...,~4(v,l))dv.       [78] 

)   tl 

Differentiating  both  sides  of  the  Eq.  [73]  v:ith  respect  to 

Mi(tl),  one  obtains 

dMj(t)     dMj(tl)    r   t  df  j(Hl(v,l) M4(v,l)) 

= +  \    dv 

dHi(tl)    dMi(tl)   J  tl         dMi(tl) 

d TT  j  ( 1 1  )    ,   t   4    df  j      d!7n  (  v  ) 

= +  \     s' dv.      [79] 

dMi(tl)    J  tl  n    diln(v)   dMi(tl) 

Because  of  the  initial  condition,  Eq.  [73],  one  has 

dM  j  (  1 1  ) 

=  dij,  1  £    j  <  3; 

dMi(tl) 

dM4(tl) 

=  0  [80] 

dMi(tl) 
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where  d  i  j  is  the  Kronecker  delta.  Differentiating  both 

sides  of  Eq.  [79],  with  respect  to  t  ,  one  obtains 

d  dMj(t)     4    dlj    dMn(t) 
=  S' . 

dt  dMi(tl)    n   dMn(t)  dMi(tl) 

In  the  matrix  form,  Eq.  [81]  can  be  written  as 


[31] 


d 
dt 


d!il(t) 
dMi(tl) 

• 

d~4 ( t ) 
L  dlii(tl)  J 


df  1 
dM  1 ( t ) 

d?4 
L  dTTl(t) 


df  1 
dM 4 ( t ) 


df4 
dU4 ( t ) 


dMl(t) 
dKi(tl) 


d!!4(t) 
L  dMi(tl)  J 


82] 


According  to  Eq.  [26],  the  solution  of  Eq.  [82]  is 


dKl(tl) 

d:ii(ti) 

dK 4 ( 1 1 ) 
dMi(tl) 


di;i(tf ) 

dMi(tl) 

* 

d ~4  (  t  f  ) 
L  dlii(tl) 


(tf.tl) 


[83] 


where  the  matrix  K(tf.tl)  is  the  solution  of  the  following 
equation : 

d?l  dll 

dMl(t)       dH4(t) 


dK(t,tl) 
dt 


K(t,tl) 


[84] 


df4        df4 
L  d~l(t)       dH4(t)  J 
with  the  initial  condition  K(tl,tl)  =  I.  Because  of 
Eq.  [30],  Eq.  [83]  can  be  written  as 
dKj(tf) 


dh'i(tl) 


=  Kji(tf.tl),   1  <  j,  i  <  3; 
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dH4(tf) 

=0  [85] 

dMi(tl) 

where  Kji(tf,tl)  is  the  jth  row,  ith  column  element  of  the 

matrix  K(tf,tl).  According  to  Eq.  [29],  one  can  have 

d'!i(tl  ,l,U  +  qS)    r   1 1  3       df  n  (  v  ,  1  ,  IJ  +  qS  ) 

=  \    dv  S(v)  S'  Kin .      [86] 

dq  J   tO         n  d'J 

Substituting  Eq.  [3  5]  and  Eq.  [86]  into  E q .  [77],  one  can 

get 

dJ(U+qS) 

dq 

3     dJ     3      r   tl  3       dfn(v,l,U+qS) 

=  S*  S'  -- S'  Kji  \     dv  S(v)  S'  Kin  

1   i   dMj(tf)  i      )  tO  n  dU 

tl  dfn(l)     dJ 

dv  S(v)  S'  S'  S'  ---  S'  Kji(tf , t 1 )Kin ( t 1 , v ) 

tO         In     dU    j   dMj  i 

t 

dv  S(v)G(v,U+qS)  [87] 

tO 

where 

E   3   dfn(l)  3    dJ    3 

G(v,U+qS)  =  S'  S'  S'  -----  S'  K ji ( tf  ,  1 1 )Kin( 1 1 , v ) 

In     dU    j    dMj   i 

L   3   dfn(v,l,U+qS) 

=  S'  S'  Qn(v),  [88] 

1   n        dU 

where 

3   dJ   3 
On(v)  =  S'  ---  S*  Kji(tf ,tl)Kin(tl,v).  [89] 

J   dMj  i 

In  the  matrix  form,  Eq.  [89]  can  be  written  as 
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01  (v) 
Q2(v) 
Q3(v) 


Kll  ...  K31 


K13  ...  K33 


Kll  ...  K 1 


L  K13  . . .  K33 


dJ 

dMl 

• 

dJ 
dM3 


[90] 


or  in  the  compact  form 

0(v)  =  K+(tl,v)  Y+(tf  ,tl)Jrn  [91] 

where  K""(tl,v)  satisfies  Eq.  [55],  Multiplying  hoth  sides 
of  Eq.  [55]  by  T  "*"(  t  f  ,  1 1  )  Jm  and  according  to  the  Eq.  [91], 
one  gets 

dQ(t) 

=  -fn  +  Q(t) 

dt 

where 

r  dfl     dfl  i 


[92] 


fm  = 


dMl 


dM3 


df3      df3 
dMl      dM3  J 
0    l*dw    0 
-l*dw   0   -U(t)  .  [93] 

0    U(t)    0 

Substituting  v  =  tl  into  Eq.  [91],  one  derives  the 
terminal  condition  for  Eq.  [92], 

Q(tl)  =  K+(tf ,tl)Jm.  [94] 

Because  Eq .  [87]  has  exactly  the  same  form  as 
Eq.  [30],  the  conjugate  gradient  method  will  minimize  the 
objective  function,  Eq.  [76],  as  it  does  for  Eq.  [14]. 
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The  optimization  procedure  will  be  the  same  as  that  for 
designing  a  180  degree  selective  rf  pulse.  The  only 
difference  is  the  tormina!  condition,  E  q .  [94],  which  is 
different  from  Eq.  [57]. 

To  find  the  terminal  condition,  E  q  .  [94],  one  can 
integrate  E q .  [84].  Because  of  Eq.  [72],  the  coefficient 
matrix  of  Eq.  [84]  will  be 
dll         dll 

—  — ™  —  —  —   #  m    9        —  — —  —  —  — 

dMl(t)      dK4(t) 


df4 
dH4(t)  ■ 


0    -0.5*l*dw    0    0 
0.5*l*dw    0      0    0 

0         0 
A(t)       B(t) 


0    0 
0    0 


[95] 


df4 
L  dMl(t) 
where 

A(t)=2*[Ml(t,l)**2  +  M2(t,l)**2  -  D(l*dw)**2]*2*Ml(t,l)**2 
B(t)=2*[Ml(t,l)**2  +  M2(t,l)**2  -  D(  l*dw)**2  ]*2*M2(  t ,  1  )**2  . 

[96] 
The  initial  condition  for  Eq.  [  ?  4  ]  is 
dHl(tl)      dHl(tl) 


dMl(tl) 

d~4 ( t 1 ) 
LdMl(tl) 


dM3(tl) 

d~4  (  1 1 ) 
dM3(tl) 


1  0  0 

0  1  0 

0  0  1 

0  0  0 


[97] 


One  can  notice  that  there  is  no  conflict  between  Eq .  [97] 
and  K ( 1 1  , 1 1 )  =  I  because  there  is  no  need  for  the  fourth 
column  of  Eq.  [97]  if  Eq.  [88]  is  used  to  obtain  the 
gradient  of  the  objective  function,  Eq.     [88],  Integrating 
Eq.  [84]  and  substituting  following  Eq.  [98]  into  Eq.  [94], 
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Jc  = 


Rl*A(tf ) 
Rl*B(tf) 

0 

1 

one  can  net  the  terminal  condition  for  0  ( t )  of  P. q .  [92], 
The  gradient  of  the  objective  function  can  he  found  by 
E q .  [99]  as  in  the  f o 1 1 o w i n g : 


[98] 


G(t.U)  =  S'  (-i!3(tfl)02(t,l)  +  M2(t,l)Q(t,l)). 


[99] 


In  summary,  the  procedure  for  designing  the  9(5  degree 
selective  rf  pulses  is  in  the  following: 

1 .  Guessing  an  initial  90  degree  rf  pulse,  for  example, 
a  90  degree  truncated  sine  function  rf  pulse,  integrate 

•7  q  .  [66]  under  the  initial  condition,  E  q  .  [67],  and  then 
integrate  E  q  .  [72]  under  the  initial  condition,  Eg.  [73], 
and  find  the  first  value  of  the  objective  function,  J(0), 
by  Eq.  [71]. 

2.  Subtituting  Eq .  [9  5]  whose  terns,  A(t)  and  B(t),  are 
determined  by  the  solution  of  Eq.  [72]  into  E q .  [84]  and 
solving  Eq.  [84]  under  the  initial  condition,  Eq.  [97], 
one  can  find  the  terminal  condition  for  Eq.  [92]  by 

Eq.  [94],  Integrate  Eq.  [92]  backwards  in  time  from  tl  to 
tO  to  get  the  vector  function  Q(t). 

3 .  Construct  the  gradient  of  the  objective  function  by 
Eq.  [99]  using  the  solution  of  Eq.  [66]. 

4.  Using  Eq.  [32]  and  Fq.  [33],  construct  a  new  rf  pulse, 
U(t).  The  numerical  coefficient,  q,  is  chosen  to  minimize 


68 


the  objective  function  J  for  given  U ( i  f  t )  and  S(i,t)  and 
the  minimizing  q  can  be  found  by  usinj;  the  Golden  Section 
Search.  The  new  rf  pulse  is  used  to  find  a  new  value  of 
the  objective  function,  J(n).  Comparing  the  new  value  with 
the  old  value  of  the  objective  function,  J(n-l),  stop  if 
J(n)  >  J(n-l),  or  construct  a  new  gradient  and  continue  by 
using  Eq.  [34]  if  J(n)  <  J(n-l). 

The  four  steps  above  can  be  inplenented  by  a  !7ortran 
program  that  is  given  in  the  Appendix  H.  The  JOB  program 
of  Appendices  B  and  H  is  given  in  the  Appendix  I.  The 
parameter  values  used  in  the  pro gran  of  Appendix  H  to 
obtain  the  90  degree  optimal  rf  pulses  in  Fig.  20,  Fig.  25 
and  Fig.  2 S  of  Chapter  VI  are  given  in  a),  b)  and  c)  of 
Appendix  J . 


CHAPTER  VI 
EXPERIMENTAL  STUDY  OF  THE  OPTIMAL  90  DEGREE 
SELECTIVE  RF  PULSES 

Before  starting  the  computation,  one  should  choose  a 
desired  transverse  magnetization  and  an  initial  trial 
function.  A  three  zero  truncated  sine  function,  as  shown 
in  Fig.  20  (dotted  line),  is  chosen  as  an  initial  trial 
function  for  the  conjugate  gradient  algorithm.  The 
response  to  the  90  degree  three  zero  sine  function  rf 
pulse  is  shown  in  Fig.  21  (dotted  line).  If  the  slice 
thickness  of  a  desired  response  is  identical  with  that  of 
the  response  to  the  three  zero  sine  function  rf  pulse  and 
if  one  chooses  WTR  =  1/6,  the  corresponding  optimal  rf 
pulse  is  shown  in  Fig.  20  (solid  line)  and  the  response  to 
the  optimal  rf  pulse  is  shown  in  Fig.  21  (solid  line). 
From  the  computer  simulation  results  one  can  conclude  that 
the  optimal  selective  90  degree  rf  pulse  has  better 
selectivity  than  the  sine  function  selective  90  degree  rf 
pulse.  There  are  less  artifacts  on  the  slice  top  and  less 
average  artifacts  outside  of  the  slice. 

The  above  conclusion  from  the  computer  simulations 
has  been  proved  by  the  experiments.  The  experiments  are 
similar  with  that  of  the  optimal  130  degree  selective  rf 
pulses.  The  instrument  used  in  the  experiment  is  the  GE 
CSI-II  2T  spectrometer/imager  system.  The  sample  is  a  tube 

69 


70 


epn)!|duiv    •a|IB|»h 


CO          U-i 

0)  o 

CO  <4-i    u 

CO  O 

cu  o 

CU 

rH 

c       e 

w 

3    0) 

^    4J     O 

rH 

C   > 

U     CO    -H 

3 

•H 

•nr  4J 

a. 

<4-l    4-1 

x:  4-i  u 

u   u 

4J             C 

<4H 

V 

-C    3 

l-l 

C  i-l 

0)     4J    Vl_ 

o  cu 

O   -H 

I— 1 

•H    CQ 

•H    >    O 

CO 

4J 

i-4             C 

E 

O  i-i 

CO    rH    -H 

•H 

C    CO 

CO     CO 

4-1 

■— s 

3    E 

•a  o 

C- 

to 

<4-l   -H 

cy  -h  o 

o 

E 

4-1 

U     4-'      U 

O    Cl 

•H    C    0) 

cu 

C    O 

CO    CU    N 

sz 

e 
E 

•H 

OJ  -o 

4-1 

CO     C0T3    -H     0) 

c 

CL) 

c 

P 

O  -H 

O     CO     1_ 

•H 

Wi  -o 

£-HXi 

co 

o  c 

H           4-» 

4-1 

N    O 

(11 

X) 

a. 

•   CO    cu 

o 

CU    CO 

*-^  rH     0) 

CU    01 

0)    3    l-i 

o 

>-.    u 

c  o.  eo  4-i 

J=     U 

•H          cu 

4J     O 

rH  HH  ~a 

« 

o 

u 

H 

CU 

-a       c 

"y 

0)    c 

•HHON 

i-  JC 

rH     CO 

Ih 

CC   4-1 

O     E     0) 

o 

cu 

n  -h  r 

4-1 

•a  -o 

v-'  4J    4J 

(1) 

c 

C 

E 

o  co 

CU    o    o 

co 

ON 

CO            4-1 

(-. 

^~v 

rH     CU 

CO 

<    CO 

3  x:  cu 

Q. 

e 

a,  4J  co 

•  -H 

c 

CU 

O  rH 

<4H      C     O 

x: 

CN 

U    iH     CH 

■a 

co   co 

• 

•   <D 

a  4->   cu 

•  vD 

00  4-1 

a;  XI    U 

CU  "»«. 

•H    4-> 

u  o 

co  >—i 

b    O 

CO        cu 

T-\ 

•o 

tl    0£ 

3     CO 

N-^ 

no    4J   4-1 

D.-H 

71 


-  o 


co 

<D 

o 

co 

4-) 

c 

o 

co 

O.XJ 

co 

c 

id 

o 

u 

c 

CO 

V 

<u 

JC 

u 

4-> 

o 

<4-l 

o 

o 

Q) 

tn 

c 

4-1 

•H 

l-l 

rH 

>N 

3 

•o 

o 

c 

9 

0) 

•H 
.-H 

3 

O 

a 

C 

CO 

• 

o 

w 

•H 

<u 

IL. 

4J 

X! 

co 

H 

i-H 

3 

• 

E 

O 

•H 

CN        • 

CO 

o 

•  CN 

M 

00 

cu 

•H        • 

4J 

U*      00 

3 

•H 

3-y-i  Cxh 

E 

o 

O 

K-l 

u 

W     O 
CU 

a) 

CO     0) 

£ 

■H      C 

H 

3    -H 

f-H 

y-<  T3 

CN 

^    -H 
i-H 

• 

0)    o 

oojr    co 

•H 

4J 

&H 

O  J5 

4-1     4J 

72 


of  water.  The  tube  length  is  100  mm  and  the  tube  diameter 
is  10  mm.  A  spin  echo  pulse  sequence,  as  shown  in  Fig.  22 
is  used  in  the  experiment.  In  the  pulse  sequence,  the  90 
degree  rf  pulse  is  frequency  selective  and  the  ISO  degree 
refocussing  rf  pulse  is  not  selective.  In  the  experiment, 
the  non-selective  ICO  degree  pulse  length  is  330  us,  the 
three  zero  sine  function  90  degree  selective  rf  pulse 
length  is  2  ms,  the  five  zero  sine  function  selective  90 
degree  rf  pulse  length  is  2  ms  as  well.  The  selective  90 
degree  rf  pulse  lengths  can  be  reduced  by  increasing  the 
output  power  level  of  the  rf  pulse  power  amplifier.  The 
echo  time,  TE,  is  30  ms  and  the  acquisition  data  block 
size  is  1024  points.  If  the  sample  axis  is  along  the  main 
magnetic  field  (z-direction),  the  chosen  project  direction 
is  the  z-direction.  One  can  make  the  experimental  rf 
pulses  by  inputing,  point  by  point,  the  numerical  data 
corresponding  to  the  rf  pulse  shapes  in  Fig.  20.  The 
experimental  results  are  shown  in  Fig.  23  and  Fig.  24. 
Comparing  Fig.  23  and  Fig.  24  with  Fig.  21,  one  can  see 
that  the  experimental  results  are  very  close  to  the 
computer  simulation  results. 

If  the  thickness  of  a  desired  slice  is  identical  with 
that  of  the  response  to  the  three  zero  sine  function  rf 
pulse  and  if  one  reduces  WTR,  that  is,  one  tries  to  get 
better  slice  quality,  for  example,  one  chooses  WTR  =  1/10, 
the  optimal  rf  pulse  obtained  is  shown  in  Fig.  25  (dotted 
line).  The  corresponding  response  to  the  dotted  line  rf 
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pulse  of  Fig.  25  is  shown  in  Fig.  26  (dotted  line).  One 
can  notice  thnt  even  though  there  are  sharper  edges  on  the 
slice,  there  are  more  artifacts  on  the  top  of  the  slice 
and  outside  of  the  slice.  The  above  conclusion  has  been 
proved  by  the  experimental  result  as  shown  in  Fig.  27. 
Thus,  one  nay  conclude  that  for  a  particular  pulse  length 
there  is  a  minimum  value  of  WTR ,  represented  by  WTRmin,  at 
a  certain  slice  thickness.  If  WTR  >  WTRmin,  one  can  always 
find  an  rf  pulse  to  which  the  response  will  be  very  close 
to  the  desired  slice.  If  WTR  <  WTRmin,  one  can  get  a 
sharper  slice  edge,  but  the  price  is  that  there  are  nore 
artifacts  on  the  top  of  the  slice  and  outside  of  the  slice. 
From  the  discussion  above  one  may  conclude  that  one  cannot 
make  a  transverse  magnetization  slice  as  sharp  as  one 
wants  with  a  single  selective  rf  pulse  at  a  certain  slice 
thickness.  All  above  conclusions  are  similar  with  that  of 
the  130  degree  selective  rf  pulses  in  the  Chapter  IV. 

If  one  tries  to  improve  the  slice  quality  for  a 
particular  pulse  length,  that  is,  to  reduce  the  value  of 
WTR,  one  must  choose  a  desired  slice  that  has  a  broader 
frequency  range.  For  example,  one  can  choose  a  desired 
slice  that  has  the  slice  thickness  which  is  identical  with 
that  of  the  response  to  a  five  zero  sine  funciton  rf  pulse 
and  WTR  =  1/10,  the  corresponding  optimal  rf  pulse  is 
shown  in  Fig.  28  (solid  line).  The  response  to  the  optimal 
rf  pulse  is  shown  in  Fig.  29.  The  experimental  results  are 
shown  in  Fig.  30  and  Fig.  31. 
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Fron  the  computer  simulation  results  and  the 
experimental  results  one  may  notice  that  the  improvements 
for  the  selectivity  of  the  90  degree  selective  rf  pulse 
after  the  optimization  procedure  are  not  so  obvious  as  for 
that  of  the  180  degree  selective  rf  pulse.  The  reason  may 
he  attributed  to  the  nature  of  the  magnetic  resonance  spin 
system,  hut  not  the  optimization  procedure. 

For  the  convenience  for  some  readers,  portions  of  the 
Fourier  coefficients  of  the  Fourier  series  for  the  90 
degree  optimal  rf  pulses  in  Fig.  20  and  Fig.  28  are  given 
in  the  Appendices  K  and  L.  The  first  20  coefficients  of 
both  cosine  and  sine  columns  in  the  Appendices  K  and  L 
will  reconstruct  the  rf  pulses  in  Fig.  20  and  Fig.  2S, 
respectively . 

The  discussion  of  the  influence  of  the  gradient 
magnitude  on  the  slice  thickness  is  similar  as  that  of 
Chapter  IV.  The  pulse  sequence  used  in  the  experiment  is 
shown  in  Fig.  22.  The  sample  is  a  tube  of  water  with  a  1 
mm  thick  bloc!;  in  the  middle  of  the  tube.  The  experimental 
result  is  shown  in  Fig.  32.  One  can  see  that  Fig.  18  and 
Fig.  32  are  quite  similar;  thus,  the  conclusions  from 
Fig.  18  and  Fig.  32  should  be  the  same  and  it  is  already 
presented  in  Chapter  IV. 
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CHAPTER  VII 
CONCLUSION 

There  are  more  than  30  articles  devoting  to  the 

selective  rf  pulse  design  problem,  except  the  articles 

which  utilize  the  optimal  control  method,  since  the 

existence  of  the  magnetic  resonance  imaging.  However, 

none  of  them  completely  solved  this  problem.  From  the 

discussion  above,  one  can  see  that  the  conjugate  gradient 

method  has  completely  solved  the  problem  of  designing  a 

selective  rf  pulse.  The  conjugate  gradient  method  is  a 

established  method  and  has  been  utilized  to  solve  the 

optimal  control  problem  since  1967  (10).  Unfortunately, 

this  powerful  method  is  just  recently  introduced  into  MRI 

and  used  to  solve  the  pulse  shape  design.  This  fact  has 

demonstrated  the  importance  again  that  the  researchers  in 

different  fields  should  understand  each  other  and  learn 

each  other. 
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APPENDIX  A 
FOURTH  ORDER  RUNGE-KUTTA  METHOD 

Because  the  fourth  order  Runge-Kutta  method  in  a 

standard  textbook,  for  example  in  (17),  is  hard  to  nrogram 

for  our  purpose,  we  rearrange  it  as  in  the  following:  if  a 

differential  equation  is 

x'(t)  =  fl(t,x,y,z,w) 

y'(t)  =  f2(t,x,y,z,w) 

[  Al  ] 
z'(t)  =  f3(t,x,y,z,w) 

w'(t)  =  f4(t,x,y,z,w) 

where  the  prime,  ',  means  differentiate  with  respect  to 

the  time  variable,  t,  the  equivalent  difference  equation 

to  above  equation  for  the  numerical  iteration  calculation 

is 

x(n+l)  =  x(n)+(Al+2*Bl+2*Cl+Dl)/6 

y(n+l)  =  y(n)+(A2+2*B2+2*C2+D2)/6 

[  A2  ] 
z(n+l)  =  z(n)+(A3+2*B3+2*C3+D3)/6 

w(n+l)  =  w(n)+(AA+2*B4+2*C4+D4)/6 

where  the  asterisk,  *,  means  multiplication  and 

Ai=dt*fi(t(n),x(n),y(n),w(n)) 

Bi=dt*fi(t(n)+dt/2,x(n)+Al/2,y(n)+A2/2,z(n)+A3/2,w(n)+A4/2) 

Ci=dt*fi(t(n)+dt/2,x(n)+Bl/2,y(n)+B2/2,z(n)+B3/2,w(n)+B4/2) 

Di=dt*fi(t(n)+dt,x(n)+Cl,y(n)+C2,z(n)+C3,w(n)+C4) 

where  dt  is  the  time  increment  and  i  =  1,  2,  3  and  4. 
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FORTRAN  PROGRAI 


APPENDIX  B 
FOR  180  DEGREE 


PULSE  DESIGN 


C  SEPT  29  19 
REAL  J 
DIMENS 
DIMENS 

DIMENS 
DIMENS 
DIMENS 
DIMENS 
DIMENS 
DIMENS 
DIMENS 


DD=f 
ML-f 
D=bl 
11,1 
DRTA 
RFA  = 
COE= 
COEF 
IDD= 
IDE= 
C  IDF= 
C  IDG  = 
C  CD  f 


C 

c 
c 
c 
c 
c 
c 
c 
c 
c 


req  .  l 
req  .  s 
ank;  P 
2,13, G 

,DRTA1 
sectio 
propor 
=renew 
ID.  of 
1  fron 

1  for 

2  no  o 
orces 
READ(1 
READ(1 
READ(1 
READ(1 
RAED(1 
READ(1 
READ(1 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
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l.J 

ION 
ION 
ION 
ION 
ION 
ION 
ION 
ION 
ION 
ION 
ncr 
tep 
=  pa 
I  d 
=  co 
n  1 
tio 
ed 
de 
si 
180 
pti 
pro 
1-* 
1 
1 
1 
1 


2,JJ,J,HUA,JA,JB,JC,JD 
X(400),Y(400),Z(400), 

RF1A(400),RF1R(400),R 

RF2A(400),RF2B(400),R 

RF0A(400),RF0B(400) ,R 

XX(400),YY(400),ZZ(40 

0MAA(400),0MBB(400),0 

W(400) ,UW(400) ,TW(200 

X1(200),Y1(200),Z1(20 

FILOUl(1200),FILOU2(l 

FIL0U4(800),FIL0U5(12 

enent;  DT=tirae  increme 

No.;  0M=func.  amplitu 

rameter  of  sine  func.  ; 

etermine  desired  Mz; 

nvergence  tolarence , DR 

ength  of  golden  sectio 

n  factor  for  integral, 

section  length,  C0EF=0 

sired  prof ile , (  1  ,  2  )  ; 

nc  func.,IDE=2  from  sa 

degree,  IDF=3  for  90 
mization;  KD=12~20,I4= 


,JQ,L1,L2,L3 
S1(400),FILI 
F1C(400),DES 
F2C(400),DES 
F0C(400),S(4 
0),G(400),GG 
MCC(400) 
),TZ(200),Z2 
0),X2(200),Y 
200),FIL0U3( 
00),FIL0U6(1 
nt;  H=initia 
de;  M=ID.  of 
NL=time  ste 


,L4,II 

NP(400) 

FIL(SOO) 

(200) 

00) 

(400) 

(200) 
2(200) 
1200) 
200) 
1  Mz; 

func . ; 
p  No.; 


run  KD  times,  each 


TA=10  50,DRTA1=1  10; 
n,RFA=0.05"0.1E-03; 

C0E=0.l"l.0E-05; 
.l'l.OE-03; 

ved  func . ; 

degree ; 

5~10, 

time  14  iterations; 


DD,DT,H,ML 

0M,M,D,P,NL 

II, 12, 13, 14, GI 

DRTA , RFA , C0E , COEF , DR 

IDD, IDE, IDF, IDG 

RX1,RX2,RY1,RY2,RZ1, 

KC, RATIO, KD 
DD,DT,H,ML 
0M,M,D,P,NL 
II, 12, 13, 14, GI 
DRTA , RFA , C0E , COEF , DRTA 
IDD, IDE, IDF, IDG 
RX1,RX2,RY1,RY2,RZ1,RZ 
KC, RATIO, KD 


TA1 
RZ2 


C  rf  p 


N3=NL-1 

M2=3*N3 

ulse  from  saved  func 

IF  (IDE.EQ.2)  THEN 
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10 


REAP (10,*)  (IDUM,FILINP(L),L=1 

DO  10  1=1, N3 

RFOA(I)«FILINP(I) 

RF0B(I)=FILINP(I+N3) 

RF0C(I)-FILINP(I+2*N3) 

CONTINUE 

ENDIF 


M2) 


Nl=(NL-l)/2+l 
N2-N1-1 

na=:u-2 

RH1-N1 

Ml=(KL+l)/2 

RM1-M1 

desired  Hz 

DO  11  1-1,11 

DESFIL(I)=H 

11  CONTINUE 

IF  (IDD.EQ.l)  THEN 
DESFIL(I1+1)=0.0 
DESFIL(Il+I2+2)=0.0 
DO  12  I-.I2 
DESFIL(I+I1+1)=- CI 

12  CONTINUE 

DO  13  1=1,13 
DESFIL(I+I1+I2+2)=H 

13  CONTINUE 
ENDIF 


16 


17 


IF  (IPD.EQ.2)  THEN 

DESFIL(Il+l)=GI/2 

PESFIL(Il+2)-0.0 

DESFIL(Il+3)»- GI/2 

DESFIL(Il+3+I2+l)=- GI/2 

DESFIL(Il+3+I2+2)=0.0 

D3SFIL( I 1+3+1 2+3) -GI/2 

DO  16  1=1,12 

DESFIL(Il+3+I)=-GI 

CONTINUE 

DO  17  1-1,13 

DESFIL(Il+3+I)=H 

CONTINUE 

ENDIF 


DO  9  I=ML+1,4*ML 
DESFIL(I)=0.0 
9  CONTINUE 


WRITE(18, 1200)  (L.DESFIL(L) ,L=1 ,4*HL) 
1200  F0RHAT(I5,E14.6) 

DO  14  1=1, ML 

DES(I)=DESFIL(I) 
14  CONTINUE 
!  I X  =  I D  .  of  gradient  change  times,  related  to  KD 
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IK=1 

112  BETA-0.0 
C  K = I D .  of  one-dimen.  optini.  iter,  times,  rela.  to  14 

K=l 

DO  15  N=1,NL 

S(N)=0.0 

S1(N)=0.0 
15  CONTINUE 
C 

111  DO  100  L=1,ML 

RL=L 

D0M=(RL-RM1)*DD 

X(1)=0.0 

Y(1)=0.0 

Z(1)=H 

U(1)=0.0 

c 

DO  101  N-1.N3 

EN=N 

IF  (K.EQ.l)  THEN 
C 

IF  (IDE.EQ.l)  THEN 
C  sine  func. 

IF  (rI.EQ.3)  THEN 

IF  (N.E0.N1)  THEN 

OMA-OM 

ELSE 

EXA=P*(RN-RN1)*DT 

OMA-OM*SIN ( EX A ) / ( EXA ) 

ENDIF 

EXB=P*((RN-RNl)*DT+DT/2) 

0MB=0M*SIN(EXB) / ( EXB ) 

IF  (N.EQ.N2)  THEN 

OMC-OM 

ELSE 

EXC=P*((RN-RN1)*DT+DT) 

OHC-OM*SIN ( EXC ) / ( EXC ) 

ENDIF 
C 

ENDIF 
C 

ELSEIF  (IDE.E0.2)  THEN 
C  a  saved  func.  as  initial  trial  func. 

0MA=RF0A(N) 

0MB=RF0MN) 

OMC-RFOC(N) 

ENDIF 
C 

OMAA(K)-OMA 

OMBB(N)-OMB 

0:iCC(N)=0riC 
C 

ELSELF  (K.NE.l)  THEN 

0MA»0MAA(N) 
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OMB=OMBB(N) 
OMC=OMCC(N) 
ENDIF 

IF  (K.EQ.l)  THEN 

RF1A(N)-0MA 

RF1B(N)=0MB 

KFlC(N)=0'iC 

ENDIF 

A1=DT*D0M*Y(K) 

B1=DT*(-D0M*X(N)-0MA*Z(N)) 

C1=DT*0MA*Y(N) 

D1=DT*(Z(N)-DES(L))**2 

A2=DT*D0M*(Y(N)+Bl/2) 

B2=DT* ( -DOM* ( X ( N ) +A 1 / 2 ) -OMB* (Z(N )+Cl /2 ) ) 

C2=DT*0MB*(Y(N)+Bl/2) 

D2=DT*(Z(N)+Cl/2-DES(L))**2 

A3=DT*D0H*(Y(N)+B2/2) 

B3=DT*(-D0M*(X(N)+A2/2)-0ME*(Z(N)+C2/2)) 

C3=DT*0MB*(Y(N)+B2/2) 

D3=DT*(Z(N)+C2/2-DES(L))**2 

A4=DT*D0M*(Y(N)+B3) 

B4=DT*  ( -DOM*  (  X  (  N  )  +A3  )  -0!'C*  (  Z  (  N )  +  C3  )  ) 

C4=DT*0MC*(Y(N)+B3) 

D4=DT*(Z(N)+C3-DES(L))**2 

X(N+l)=X(N.)  +  (Al+2*A2+2*A3+A4)/6 
Y(N+l)=Y(N)+(Bl+2*B2+2*B3+B4)/6 

Z(N+l)=Z(N)+(Cl+2*C2+2*C3+C4)/6 

W(N+l)=W(N)+(Dl+2*D2+2*D3+D4)/6 
101  CONTINUE 

IF  (K.EQ.l)  THEN 

IF  (L.EQ.H1)  THEN 

PRINT  *,  Z(l) ,Z(N1),Z(NL) 

ENDIF 

ENDIF 

IF  (X.EO.l)  THEN 

X1(L)=X(NL) 

Y1(L)=Y(NL) 

Z1(L)=Z(NL) 

ENDIF 

TZ(L)«Z(NL) 

TW(L)=U(NL) 
100  CONTINUE 

J0=0.0 
JA=0.0 

DO  21  L=1,ML 

J=RZ2*(TZ(L)-DES(L))**2+RZl*TTv.'(L) 
JO=JO+J 

J*RZ2*(TZ(L)-DES(L))**2 
JA=JA+J 
21  CONTINUE 
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C 

IF  (K.EQ.l)  THEN 

PRINT  *,  JO, J A 

ENDIF 
C  no  optimi.  goto  the  end 

IF  (IDG.E0.2)  THEN 

GOTO  402 

ENDIF 
C 

DO  200  L=1,ML 

RL=L 

do;:=(rl-r;h)*dd 

X(1)=0.0 

Y(1)=0.0 

Z(l)-H 

W(1)=0.0 

c 

DO  201  N=1,N3 

RN-R 

IF  (K.EQ.l)  THEN 


C 

C 


c 

c 
c 


IF  (IDE.EQ.l)  THEN 

IF  (M.EQ.3)  THEN 

IF  (N.EQ.N1)  THEN 

OMA-OM 

ELSE 

EXA=P*(RN-RN1)*DT 

OMA=OM*SIN(EXA) / ( EXA ) 

ENDIF 

EXB=P*((RN-RNl)*DT+DT/2) 

OMB=OK*SIK(EXB)/(EXB) 

IF  (N.EQ.N2)  THEN 

OMC=OM 

ELSE 

EXC=P*((RN-RN1)*DT+DT) 

OMC=OM*SIN ( EXC ) / ( EXC ) 

ENDIF 

ENDIF 

ELSEIF  (IDE.E0.2)  THEN 

OMA=RFOA(N) 

OM'B=RFOB(N) 
OMC=RFOC(N) 

ENDIF 

OMAA(N)-OMA 
OMBB(N)-OMB 

OMCC(N)=OMC 

ELSELF  (K.NE.l)  THEN 
OMA=OMAA(N) 


OMB=OMBB(N) 
OMC=OMCC(N) 

ENDIF 


93 


IF  (K.EQ.l) 
RF1A(N)=0MA 
RF1B(N)=0MB 
RF1C(N)=0MC 
ENDIF 


THEN 


AUDT' 
E1=DT: 
C1  =  DT- 
D1=DT* 
A2  =  DT: 
B2=DT- 
C2  =  DTJ 
D2=DTS 
A3  =  DT* 
B3=DT< 
C3  =  DT! 
D3=BTS 
A4=DT- 
B4  =  DT; 
C4=DT" 
D4=DT> 


DOH*Y(N) 

(-DOM*X(N)-OMA*Z(N)) 

OMA*Y(N) 

(Z(N)-DES(L))**2 

D0M*(Y(N)+Bl/2) 

( -DOM* ( X ( N ) +A 1 / 2 ) -OHB* ( Z ( N ) +C 1 / 2 ) ) 

OMB*(Y(N)+Bl/2) 

(Z(N)+Cl/2-DES(L))**2 

DOM*(Y(N)+B2/2) 

(-DOM*(X(N)+A2/2)-OMB*(Z(N)+C2/2)) 

OMB*(Y(N)+B2/2) 

(Z(N)+C2/2-DES(L))**2 

DOM*(Y(K)+B3) 

( -DOM* ( X ( N ) + A3 ) -OMC* ( Z ( N ) +C3 ) ) 

OMC*(Y(N)+B3) 

(Z(N)+C3-DES(L))**2 


)+(Al+2*A2+2*A3+A4)/6 

)+(Bl+2*B2+2*E3+B4)/6 
)+(Cl+2*C2+2*C3+C4)/6 
)+(Dl+2*D2+2*D3+D4)/6 


X(N+1)=X(N 
Y(N+1)=Y(N 
Z(N+1)=Z(N 
W(N+1)=W(N 
201  CONTINUE 
XX(1)=0.0 
YY(1)=0.0 

ZZ(l)  =  2*RZ2*(rfZ(L)-DES(L)  ) 
WU(1)=RZ1 


C  in 
C 


BO  202  N=lfN3 

RN  =  N 
tegrating  adjoint 

DDT=-DT 


enua.  back  wards  in  tine 


IF  (K.EQ.l)  THEN 
IF  (IDE. EG. 1)  THEN 


IF    (M. 

EQ.3)    THEN 

IF    (N. 

EQ.N1)    THEN 

o:iA=Oi- 

! 

ELSE 

EXA=P* 

!(RN-RN1)*DT 

0MA=0I 

I*SIN(EXA)/(EXA) 

ENDIF 

EXB=P* 

:((RN-RNl)*DT+DT/2) 

OMB=OM*SIN ( EXB) / ( EXB ) 
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C 
C 
C 


IF  (N.EQ.N2)  THEN 

OMC-OM 

ELSE 

EXC=P* ( ( RN-RH 1 )*DT+DT) 

OIiC=OM*SIN(EXC)/(EXC) 

ENDIF 

ENDIF 

ELSEIF  (IDE.  EC).  2)  THEN 


OMA-RFOA(N) 

OMB=RFOE(N) 

OMC-RFOC(K) 

ENDIF 
C 

ELSELF  (K.NE.l)  THEN 

OMA-OMAA(N) 

OMB-OKBB(N) 

OMC-OMCC(K) 

ENDIF 
C 

ZD-(Z(N)+Z(N+l))/2 

Al-DDT*DOM*YY(N) 

3 1=DDT* ( -DOM*  XX ( N ) -OHA*ZZ ( N ) ) 

C 1  =  DPT* ( OMA*YY ( N ) -2* ( Z ( N ) -DES ( L ) ) *WW (N ) ) 

Dl-0.0 

A2-DDT*DOM* ( YY ( N ) +E 1 / 2 ) 

B2=DDT* (-DOM* ( XX ( N ) +A 1 / 2 ) -OHB* (ZZ( N ) +C 1 / 2 ) ) 

C2=DDT*(OMB*(YY(N)+Bl/2)-2*(ZD-DES(L))*(WW(N)+Dl/2)) 

D2=0.0 

A3-DDT*DOM* ( YY(N )+B2/2) 

B3=BDT*(-DOK*(XX(N)+A2/2)-OMB*(ZZ(N)+C2/2)) 

C3=DDT*(OMB*(YY(N)+B2/2)-2*(ZD-DES(L))*(WW(N)+D2/2)) 

D3=0.0 

A4=DDT*D0M*(YY(N)+B3) 

B4-DDT* ( -DOM* ( XX( N) +A3 )-OMC*( ZZ (N )+C3 ) ) 

C4=DDT*(OMC*(YY(N)+B3)-2*(Z(N+l)-DES(L))*(WW(N)+D3)) 

D4=0.0 
C 

XX(N+l)=XX(N)+(Al+2*A2+2*A3+A4)/6 

YY(N+l)=YY(N)+(Bl+2*B2+2*B3+B4)/6 

ZZ(N+l)=ZZ(N)+(Cl+2*C2+2*C3+C4)/6 

WW(N+l)=WW(N)+(Dl+2*D2+2*D3+D4)/6 
202  CONTINUE 
C  G(N)=gradient  of  objective  func. 

DO  23  N-l.NL 

G(N)=-Z(N)*YY(NL+1-N)+Y(N)*ZZ(NL+1-N) 
2  3  CONTINUE 

IF  (L.EQ.l)  THEN 

DO  24  N=1,NL 

GG(N)=0.0 
24  CONTINUE 

ENDIF 
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DO  25  N-l.NL 
GG(N)-G6(N)+G(N) 
2  5  CONTINUE 
C 

200  CONTINUE 
C 

SC1=0.0 

DO  31  1=1, N2 

SG1=SG1+(C0E*GG(2*I-1))**2 

31  CONTINUE 
SG2=0.0 

DO  32  I-1.K4 
SG2=SG2+(C0E*GG(2*I))**2 

32  CONTINUE 
DOT2=((COE*GG(l))**2+(COE*GG(NL))**2+4*SGl+2*SG2)*DT/3 
IF  (K.NE.l)  THEN 

BETA=D0T2/D0T1 

END  IF 

D0T1=D0T2 
C 

DO  33  N=1,NL 

S ( N ) =-GG ( N)+BETA*S ( N ) 
5  3  CONTINUE 
C  initial  setting  for  golden  section  search 

HUA=0. 618034 

RFA1  =  RFA*H1JA 

RFA2=RFA 

L1=RFA*KUA**2 

L2=RFA*HUA 

L3=0.0 

L4=RFA 

RLN=L1-L3 

RLM=L4-L2 
C  after  KC  times  G(N),  setting  is  changed  to  keep  converge 

IF  (K.NE.l)  THEN 

IF  (K.GT.KC)  THEN 

C0EF=RATI0*C0EF 

ENDIF 

RFA1-C0EF*RFA*HUA 

RFA2=COEF*RFA 

L1-C0EF*RFA*HUA**2 

L2=COEF*RFA*HUA 

L3=0.0 

L4«C0EF*RFA 

RLN=L1-L3 

RLN=L4-L2 

ENDIF 
C 

ID=0 

1  =  0 
C 

333  DO  300  L=1,KL 

RL«L 

D0M=(RL-RM1)*DD 


c 
c 
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X(1)=0.0 
Y(l)=U.O 
Z(l)-H 

','(1)  =  0.0 

DO  301  K-1,N3 

RN=N 

IF  (K.EQ.l)  THEN 

IF  (IDE.EQ.l)  THEM 

IF  (K.EQ.3)  THEN 

IF  (N.EQ.N1)  THEN 

0MA=0K+RFA1*S(N) 

ELSE 

EXA=P*(RN-RN1)*DT 

OMA-OM*SIM ( EXA) / ( EXA ) +RFA1*S( N ) 

ENDIF 

EXB=P*((RN-RNl)*DT+DT/2) 

0MB=0M*SIN(EXB)/(EXB)+RFAl*(S(N)+S(N+l))/2 

IF  (N.E0.N2)  THEN 

OMOOM 

ELSE 

EXC=P* ( ( RN-RN 1 ) *DT+DT ) 

0MC=0M*SIN(EXC)/(EXC)+RFA1*S(N+1) 

ENDIF 

ENDIF 

ELSEIF  (IDE.EQ.2)  THEN 

OH A =RFO  A ( N ) +RF A 1 *  S ( N ) 

01  iB=RFOB (  N  )  +  RFA 1* ( S ( N )+S ( N+l )  )  /2 

0MC-RF0C(N)+RFA1*S(N+1) 

ENDIF 

ELSELF  (K.NE.l)  THEN 
OHA-OMAA ( N ) +RFA 1*S (N ) 
OMB-OMBB(  N )  +F.FA  1  *  (  S  (  N  )  +S  (  N  + 1  )  )  / 2 
0MC-OMCC( N ) +RF A 1*S ( N+l ) 

ENDIF 

ONAA(N)-OMA 
OMBB(K)-OMB 

OMCC(K)«OMC 
RF2A(N)=0MA 
RF2B(N)«0KB 
RF2C(N)»0MC 

ENDIF 

A1=DT*D0M*Y(N) 
B1=DT*(-D0M*X(N)'-0KA*Z(N)) 
C1=DT*0MA*Y(N) 
D1=DT*(Z(N)-DES(L))**2 


c 


c 
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A2  =  DT*D0.I*(Y(N)  +  31/2) 

B2=DT*  (-DOM*  (  X  ( K )  +  A 1  /  2  )  -03  IB*  ( Z  (  N )  +C 1  /  2  ) ) 

C2=DT*0MB* ( Y ( N ) +E 1 / 2 ) 

D2=DT* ( Z ( N ) +C 1 / 2-DES ( L ) ) ** 2 

A3=DT*D0M*(Y(N)+B2/2) 

B3=DT*  ( -DOM* (  X  (N )+A2/2 ) -OMB*  ( Z ( II )  +C2 / 2  )  ) 

C3=DT*0MB*(Y(N)+B2/2) 

D3=DT*(Z(N)+C2/2-DES(L))**2 

A4  =  DT*D0ii*(Y(N)  +  B3) 

B4=DT* (-DOM* ( X ( N ) +A3 ) -OMC* ( Z ( N ) +C3 ) ) 

C4=DT*0MC*(Y(N)+B3)" 

D4=DT*(Z(N)+C3-DES(L))**2 

X(N+l)=X(N)+(Al+2*A2+2*A3+A4)/6 
Y (N+l ) =Y ( N )+(Bl+2*B2+2*B3+B4 ) /6 
Z(N+l)=Z(N)+(Cl+2*C2+2*C3+C4)/6 
\\T(N+l)=W(N)  +  (Dl+2*D2+2*D3+D4)/6 

301  CONTINUE 

X2(1)=X(NL) 
Y2(1)=Y(NL) 
Z2(1)=Z(KL) 

TW(1)=W(NL) 
300  CONTINUE 

JJ=0.0 
J3=0.0 

DO  41  L-l.ML 

J=RZ2*(Z2(L ) -DES ( L ) ) **2+RZl *TW(L ) 
JJ=JJ+J 

J=RZ2* ( Z2 ( L ) -DES (L) )**2 
JB=JB+J 
41  CONTINUE 

1  =  1  +  1 

IF  (I.EQ.l)  THEN 

J2=JJ 

RFA1=RFA1*HUA 

GOTO  333 

ELSEIF  (I.EQ.2)  TEEN 

J1=JJ 

ENDIF 

IF  (ID.EQ.l)  THEN 

J1=JJ 

ELSEIF  (ID.E0.2)  THEN 

J  2= J  J 

ENDIF 

IF  (K.GT.l)  THEN 

IF  (Jl.GT.JC)  THEN 

ID0=1 

GOTO  400 

ELSEIF  (JB.GT.JD)  THEN 

ID0=1 

GOTO  400 
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EN  DIP 

EM  DIP 
C 

IF  (K.EQ.I4)  THEN 

IB0=1 

GOTO  A 00 

EKDIF 
C  p, olden  section  search  is  working 

IF  (J2.LT.J1)  THEN 

RFA2=RFA2-RLN 

ID=2 

L3=L1 

L1=L2 

J1  =  J2 

L2=L3+HUA*RFA2 

RLM-L4-L2 

RLN=L1-L3 

EFA1=L2 

GOTO  333 
C 

ELSE IF  (J2.GT.J1+DRTA)  THEN 

RFA2-RFA2-RLH 

ID=1 

L4=L2 

L2  =  L1 

J2=J1 

L1=L4-HUA*RFA2 

RLH=L4-L2 

RLN=L1-L3 

RFA1=L1 


GOTO  333 

ELSE 

JC=J1 

JD-JIi 

PRINT  *,  Jl.JA.JE 

EKDIF 

DO  43  N=1,NL 

S1(N)=S1(N)+RFA1*S(K) 

43 

CONTINUE 

K=K+1 

IF  (J2.LT.J1+DRTA1)  THEN 

ID0=1 

GOTO  400   ' 

ELSE  ID0=2 

GOTO  111 

ENDIF 

400 

CONTINUE 

IF  (IDE.EQ.2)  THEN 
IF  (IK.LT.KD)  THEN 
DO  401  1=1, N3 
RF0A(I)=RF2A(I) 
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i?FOB(I 
FOC(I 


)=I!F2:.(T) 
)=RF2C(I) 


401  CONTINUE 


IK=IK+1 
PRINT  *,IK 
GOTO  112 
ENDIF 
END  IF 

40  2  DO  4  4 
FIL0U1 
FIL0U1 
FIL0U1 

44  CONTIN 
DO  4  5 
FIL0U2 
F I LOU 2 
FIL0U2 

45  CONTIN 
DO  46 
FIL0U3 
FIL0U3 
FIL0U3 
SQ=X1( 
FIL0U3 

4  6  CONTIN 

DO  4  7 

FIL0U4 

FIL0U4 

FIL0U4 

S0=X2( 

F I LOU 4 
47  CONTIN 

M3=4*M 

DO  48 

FIL0U5 
4  3  CONTIN 

DO  49 

F I  LOU 6 
49  CONTIN 

WRITE( 

WRITE( 


WRITE ( 

V.T.ITE( 
WRITE ( 

WRITE ( 
1100  FORMAT 
STOP 
END 


KK=1,N3 

(KK)-RFIA 

(KJC+N3)=R 

(D'+2*N3) 

UE 

K?:=1,N3 

(KK)=RF2A 

(KK+iJ3)=S 

(K    +2*;'3) 

UE 

KK-l.ML 

(KK)=Z1(K 

(KK+HL)=X 

(KK+2*I5L) 

KK)**2+Y1 

(KK+3*ML) 

UE 

J\  K  =  i  9V1U 

(KK)=Z2(K 
(KK+ML)-X 
(KK+2*ML) 
KK)**2+Y2 
(KK+3*ML) 
UE 

T 

KK=1,M2 
(XK)-O.O 

UE 

KK-1.K3 

(KK)=0.0 

UE 

22,1100) 

13,1100) 

14,1100) 

15,1100) 

16,1100) 

17,1100) 

(I5.F14.6 


(KS) 


FIB(KK) 
=RF1C( KK) 


(KK) 

F2B(KK) 

=RF2C(KK) 


K) 

1(KK) 
=Y1(KK) 
(KK)**2 
=SQRT(SQ) 


2(KK) 
=Y2(KK) 
(KK)**2 
»S0RT(SQ) 


(L,FIL0U1(L),L=1,M2) 
(L,FIL0U2(L) ,L-1,M2) 
(L,FIL0U3(L),L=1,M3) 
(L,FIL0U4(L),L=1,M3) 

(L,FILOU5(L),L=l,M2) 
(L,FIL0U6(L),L=1,M3) 
) 


APPENDIX  C 
PARAMETER  VALUES  OF  THE  PEOAGRAM  IN  THE  APPENDIX  B 
TO  OBTAIN  THE  OPTIMAL  RF  PULSES 
OF  FIG.  2,  FIG.  11  AND  FIG.  14 

a)  The  parameter  values  for  the  optimal  rf  pulse  of  Fig.  2 


2111.0, 

1.  11111 1E-05, 

2000, 

159 

-26911.0, 

3, 

0.7, 

25300.0, 

91 

69, 

15, 

69, 

5, 

2000 

50.0, 

0.5E-03, 

0.5E-05, 

0.5E-03, 

10.0 

2, 

1, 

1, 

0 

1, 

1, 

1, 

1, 

1, 

60, 

0.7, 

20 

b)  The  parameter  values  for  the  optimal  rf  pulse  of  Fig.  11 


2111.0, 

1.11111 1E-05, 

2000, 

159 

-26911.0, 

3, 

0.7, 

25300.0, 

91 

70, 

17, 

70, 

5, 

2000 

50.0, 

0.5E-03, 

0.5E-05, 

0.5E-03, 

10.0 

2, 

1, 

1, 

0 

1, 

1, 

1, 

1, 

1, 

60, 

0.7, 

20 

c)  The  parameter  values  for  the  optimal  rf  pulse  of  Fig.  14 


2111.0, 

1.111111E-05, 

2000, 

159 

-39615.0, 

3. 

0.7, 

3S350.0, 

91 

64, 

25, 

64, 

5, 

2000 

50.0, 

0.5E-03, 

0.5E-05, 

0.5E-03, 

10.0 

2, 

1, 

1, 

0 

1, 

1, 

1, 

1, 

1, 

60, 

0.7, 

15 

100 


APPENDIX  D 

FORTRAN  PROGRAM  FOR  REFOCUSING  PULSE  STUDY 

C  JAN'.  8.  1987 
C  program  SELE53 

DIENSIOH  X ( 400 ), Y ( 400 ), Z ( 400 ),FILINP( 1200) 

DIMENSION  RF1A(400) ,RF1B(400) ,RF1C(400) 

DIMENSION  RF0A(400) ,RF0B(400) ,RF0C(400) 

DIMENSION  XO(IOOO) ,YO(1000) ,Z0(1000) 

DIMENSION  X 1 ( 1 000 ) , Yl ( 1 000 ) , Z 1 ( 1 000 ) 

DIMENSION  XA( 1000 ), YA ( 1 000 ),ZA( 1000) 

DIMENSION  XB ( 1 000 ), YB ( 1 000 ),ZB( 1000) 

DIMENSION  XYl(lOOO) ,XYA(1000) .XYB(IOOO) 

DIMENSION  XY0( 1000), XYR( 1000) 

DIMENSION  FILOU1(1200),  FIL0U2( 1200) , FIL0U3( 1C00) 

DIMENSION  FIL0U4(1000),  FIL0U5( 1 200) , FIL0U6( 1000) 
C  IDD  desired  profile  (1,2,3),  IDE  run(l)  run(2) 
C  IDH  for  different  rotation  axis 

READ (11,*)  DD,DT,H,ML,IHW 

RLAD(11,*)  OM,M,D,P,NL 

READ(11,*)  II, 12, 13, 14, GI 

READ (11,*)  IDD, IDE, IDF, IDG, IDH 

PRINT  *,  DD,DT,H,ML,IKW 

PRINT  *,  OM,M,D,P,NL 

PRINT  *,  II, 12, 13, 14, GI 

PRINT  *,  IDD, IDE, IDF, IDG, IDH 
C 

N3=NL-1 

M2=3*N3 

IF  (IDE.EQ.2)  THEN 

READ (10,*)  (IDUM,FILINP(L),L=1,M2) 

DO  10  1=1, N3 

RF0A(I)=FILINP(I) 

RF0B(I)=FILINP(I+N3) 

RF0C(I)=FILINP(I+2*N3) 
10  CONTINUE 

ENDIF 
C 

X(1)=0.0 

Y(l)-H 

Z(1)=0.0 

c 

DO  100  L=1,ML 
RL=L 

D0M=(RL-RM1)*DD . 
C 

DO  101  N=1,N3 
RN  =  N 


101 


102 


IF  (IDE.EQ.l)  THEN 

IF  (M.EQ.l)  THE:i 

OMA-OM 

OMB=OM 

OMC-OM 

ELSEIF  (M.EQ.3)  THEN 

IF  (N.EQ.N1)  THEN 

0MA=OM 

ELSE 

EXA=P*(RN-RN1)*DT 

0MA=0M*SIN ( EXA ) / ( EXA ) 

ENDIF 

EXB=P*((RN-RNl)*DT+DT/2) 

0;iB=OK*S IN (  EXB )  /  (  EXB  ) 

IF  (N.EQ.22)  THEN 

OMC-OM 

ELSE 

EXC-P*((RN-RN1)*DT+DT) 

OMC=OM*SIN(EXC)/(EXC) 

ENDIF 

ENDIF 

ELSEIF  (IDE.EQ.2)  THEN 

OMA-RFOA(N) 

OMB-RFOBNN) 

OMC-RFOC(N) 

ENDIF 

RF1A(N)-0MA 

RF1B(N)-0MB 

RF1C(N)=0MC 

IF  (IDII.EQ.l)  THEN 

A1=DT*D0M*Y(N) 

B 1  =  DT* ( -DOM*X ( N ) -0MA*Z( N ) ) 

C1=DT*0MA*Y(N) 

A2«DT*DOM*(Y(N)+Bl/2) 

B2=DT*(-DOM*(X(N)+Al/2)-OMB*(Z(N)+Cl/2)) 

C2=DT*OMB*(Y(N)+Bl/2) 

A3-DT*D0M*(Y(N)+B2/2) 

B3=DT*(-DOM*(X(K)+A2/2)-OMB*(Z(N)+C2/2)) 

C3=DT*0MB*(Y(N)+B2/2) 

AA=DT*D0M*(Y(N)+B3) 

B4-DT*(-DOM*(X(N)+A3)-OKC*(Z(N)+C3)) 

C4-DT*0MC*(Y(N)+B3) 

ELSEIF  (IDH.EQ.2)  THEN 

A1=DT*(D0M*Y(N)-0MA*Z(N)) 

B1-DT*(-D0M*X(N)) 

C1=DT*0MA*X(N) 

A2-DT*(D0M*(Y(N)+Bl/2)-0MB*(Z(N)+Cl/2)) 

B2-DT*(-DOM*(X(N)+A1*0.5)) 

C2»DT*0MB*(X(N)+Al/2) 

A3-DT*(D0M*(Y(N)+B2/2)-0MB*(Z(N)+C2*0.5)) 

B3=DT*(-D0M*(X(N)+A2*0.5)) 
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C3=DT*OMB*(X(N)+A2/2) 

A4=DT*( DOM* ( Y ( N )+B3 )-OMC* ( Z ( N ) +C3 ) ) 
B4=DT*(-D0M*(X(N)+A3)) 

C4=DT*0HC*(X(N)+A3) 

ELSEIF  (IDH.EQ.3)  THEN 

A1=DT*(-0MA*Z(N)) 

B1=DT*D0M*Z(N) 

C 1 =DT* ( OMA*X( N ) -POM* Y ( N ) ) 

A2«=DT* ( -OMB* ( Z ( N )  +C 1  / 2 )  ) 

B2=DT*D0M*(Z(N)+Cl/2) 

C  2  =  DT* ( OMB* ( X ( N ) +A 1 / 2  ) -DOM* ( Y ( N ) +E 1 / 2 ) ) 

A3=DT*(-0MB*(Z(N)+C2/2)) 

B3=DT*D0M*(Z(N)+C2/2) 

C3=DT*(0MB*(X(N)+A2/2)-D0H*(Y(N)+B2/2)) 

A4=DT*(-0MC*(Z(N)+C3)) 

B4=DT*D0M*(Z(K)+C3) 

C4=DT*  ( OMC* (  X  (N )  +A3  )  -DOM*  (  Y  ( 11 )  +B3  ) ) 

ENDIF 

X(N+l)=X(N)+(Al+2*A2+2*A3+A4)/6 

Y(N+l)=Y(N)+(Bl+2*E2+2*33+B4)/6 
Z(N+l)=Z(N)+(Cl+2*C2+2*C3+C4)/6 
101  CONTINUE 

IF  (L.EQ.H1)  THEN 

PRINT  *,  X(1),X(N1),X(L) 

PRINT  *,  Y(1),Y(N1),Y(NL) 

PRINT  *,  Z(1),Z(N1),Z(NL) 

ENDIF 

Z1(L)=Z(NL) 

xi(l)=x(i;l) 

Y1(L)=Y(NL) 

X Y 1 ( L )=SQRT ( XI ( L ) **2+Y 1 (L ) **2 ) 

100  CONTINUE 

DO  44  KK=1,N3 
FIL0U1(KK)=RF1A(KK) 

FILOU 1 ( KK+N3 ) =RF1B ( KK ) 
FIL0U1(KK+2*N3)=RF1C(KX) 
44  CONTINUE 

WRITE(22,1100)  (L,FIL0II1(L)  ,L=1,3*N3) 

WRITE(17,1100)  (L, XI (L), Lei, ML) 

VRITE(13,1100)  (L,Y1(L),L=1,ML) 

WRITE(19,1100)  (L,Z1(L),L=1,ML) 

URITE(20,1100)  (L,XY1(L) ,L-1,ML) 
1100  F0RMAT(I5,E14.6) 
STOP 
END 


APPENDIX  E 
FOURIER  COEFFICIENTS  FOR  THE  FOURIER  SERIES 
OF  THE  OPTIMAL  RF  PULSE  IN  FIG.  2 


REPRESEKTION 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 


Cosine 

0, 

.296246E+04 

0, 

.602884E+04 

0. 

.603475E+04 

0. 

,  623589E+04 

0, 

.382474E+04 

0. 

.934820E+04 

0, 

.690043E+04 

0. 

.359159E+04 

0. 

.256553E+04 

0. 

,  160030E+04 

0, 

.122225E+04 

0, 

.771072E+03 

0. 

■599388E+03 

0. 

.376579E+03 

0. 

.2S8086E+03 

0. 

.177351E+03 

0, 

.12S165E+03 

0. 

.758214E+02 

0, 

.546459E+02 

0. 

.213661E+02 

0. 

,  183438E+02 

0. 

.924536E+00 

0, 

.111633E+02 

0, 

.225646E+01 

0, 

•673807E+01 

0, 

.609533E+01 

0, 

.504489E+01 

0, 

.421767E+01 

0, 

.760243E+01 

0, 

.255727E+01 

0, 

.457754E+01 

0, 

.672835E+01 

0, 

.353689E+01 

0. 

.637534E+00 

0, 

.101881E+02 

0. 

.218592E+01 

0, 

.753347E+01 

0. 

.971475E+00 

0, 

.833246E+01 

0, 

,  126060E+00 

Sine 

O.OOOOOOE+00 

■0.374348E+03 

0.776027E+03 

■0.118075E+04 

0.173430E+04 

■0.310588E+04 

0.246532E+04 

■0.170867E+04 

0.136224E+04 

■0.104353E+04 

0.863594E+03 

■0.667074E+03 

0.543747E+03 

•0.424560E+03 

0.329281E+03 

-0.267328E+03 

0.183564E+03 

■0.157226E+03 

0.919727E+02 

■0.751788E+02 

629E+02 

716E+02 

393E+02 

127E+02 

695E+02 

642E+02 

206E+01 

300E+01 

217E+01 

141E+01 

544E+01 

982E+01 

891E+01 

0.383432E+01 


157 
919 
751 

0.329 
-0.357 

0.210 
■0.259 

0.109 
■0.118 

0.208 
■0.435 
-0.697 

0.506 
■0.269 
■0.547 

0.132 

18 
■0.101 

0.842 
•0.556 

0.392 
■0.399 

0.204 


133E+02 
952E+01 
266E+01 
170E+01 
435E+01 
690E+01 
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APPENDIX  F 
■OURIER  COEFFICIENTS  FOE  THE  FOURIER  SERIES  REPRESENTATION 

OF  THE  OPTIMAL  RF  PULSE  IN  FIG.  14 

Sine 

O.OOOOOOE+OO 

-0.412471E+03 
0.907366E+03 

-0.171388E+04 
0.125197E+04 

-0.616422E+03 
0.433608E+03 

-0.248337E+03 
0.155689E+03 

-0.923386E+02 
0.571649E+02 

-0.273676E+02 
0.2510R7E+02 

-0.237399E+01 
0.108717E+02 
0.189106E+01 
0.865308E+01 
0.412969E+01 
0.352519E+01 
0.347007E+01 
0.824660E+01 
0.281181E+00 
0.252998E+01 
O.S98227E+01 

-0.1S4282E+01 
0.509570E+01 
0.202300E+01 
0.351510E+01 
0.183202E+01 
0.868455E+00 
0.426226E+01 
0.137229E+01 

-0.199483E+00 
0.457037E+01 

-0.280655E+00 
0.121986E+01 
0.143490E+01 
0.135804E+01 
0.359790E+00 
0.413855E+00 
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Cosine 

1 

-0. 

291488E+04 

2 

0. 

596038E+04 

3 

-0. 

647590E+04 

4 

0. 

803627E+04 

5 

-0. 

456489E+04 

6 

0. 

165962E+04 

7 

-0. 

997046E+03 

8 

0. 

451111E+03 

9 

-0. 

25991cE+03 

10 

0. 

118195E+03 

11 

-0. 

748929E+02 

12 

0. 

23O019E+02 

13 

-0. 

275648E+02 

14 

-0. 

214526E+01 

15 

-0. 

106583E+02 

16 

-0. 

416702E+01 

17 

-0. 

679630E+01 

18 

-0. 

417355E+01 

19 

-0. 

331207E+01 

20 

-0. 

269037E+01 

21 

-0. 

352532E+C1 

22 

-0. 

15184°E+01 

23 

-0. 

142062E+01 

24 

-0 

143621E+01 

25 

-0. 

130772E+01 

26 

-0. 

287690E+00 

27 

-0. 

962  701E+00 

28 

0. 

378960E-01 

29 

0. 

641276E-01 

30 

-0. 

882932E+00 

31 

0. 

143774E+01 

32 

0. 

274126E+00 

33 

-0 

12476SE+01 

34 

0. 

308371E+01 

35 

-0 

692546E+00 

36 

0. 

518338E+00 

37 

0 

109133E+01 

38 

0 

121201E+01 

39 

0 

676590E+00 

40 

-0 

.377214E+00 

APPENDIX  G 
FORTRAN  PROGRAM  FOR  FOURIER  COEFFICIENTS 

C  NOV.  23RD  1986 

DIMENSION  A(400),B(400),F(606),FILINP(606),FS(600) 

DIMENSION  SG1(400) ,SG3(400) , SG4(400) , SG2(400) ,S(400) 

READ(11,*)  DD,DT,H,ML,ID 

READ ( 11,*)  OH , M , D , P , NL , N J 

PRINT  *,  DD,DT,H,ML,ID 

PRINT  *,  OM,M,D,P,NL,NJ 

N3-NL-1 

READ (10,*)  (IDUM,FILINP(L),L=1,N3) 

DO  10  1=1, N3 

F(I)=FILINP(I) 

10  CONTINUE 
TPI=6. 2831853 
A(1)=0.0 
B(1)=0.0 
S(1)=0.0 

DO  12  N=1,NL 

DO  11  K=1,N3 

A ( N )=A( N )+F( K ) *COS (TPI* ( N-l )*K/N3) 

B(N)=B(N)-F(K)*SIN(TPI*(N-1)*K/N3) 

11  CONTINUE 
A(N)«2*A(N)/N3 
E(N)=2*B(N)/N3 

12  CONTINUE 
A(l)=A(l)/2 

C 

DO  22  K=1,N3 
DO  21  N=1,NJ 

SS=B(N)*SIN(TP1*(N-1)*K/N3) 
S(K)=S(K)+A(N)*C0S(TPI*(N-1**K/N3)-SS 

21  CONTINUE 
FS(K)=S(K) 
FS(K+N3)=S(K) 
FS(K+2*N3)=S(K) 

22  CONTINUE 

WRITE(14,1100)  (L,A(L),L=1,NL) 

WRITE(15,1100)  (L,B(L),L=1,NL) 

WRITE ( 1 6 , 1 1 00 )  ( L , FS ( L ) , L= 1 , 3*N3 ) 

1100  F0RMAT(I5,E14.6) 

STOP 

END 
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APPENDIX  II 
FORTRAN  PROGRAM  FOR  90  DEGREE  PULSE  DESIGN 

C  Sept.  29  1986 
C  program  SELE52 

REAL  Jl,  J2,  J,  JJ.HUA.II,  JA,  JE,  JC , JD , JO , LI , L2 , L3  ,  L4 


C 
C 

C 

c 


IDD 

IDF 
IDO 


DIMEN 
DIMEN 
DIMEN 
DIMEN 

DIMEN 
DIMEN 
DIMEN 
DIMEN 
DIMEN 
DIMEN 
DIMEN 
DIMEN 

desir 
180  d 
Dz(l) 
EEAD( 
READ( 
READ( 
READ( 
READ( 


SIGN  X(201),Y(201),Z(201),S1(201),FILINP(606) 

SION  XA(201),YA(201),ZA(201),U3(160),Z2(160) 

SION  XB(201),YB(201),ZB(201),WB(201) 

SION  RF1A(201),RF1B(201),RF1C(201),DESFIL(640) 

SION  RF2A(201),RF2B(201),RF2C(201),DXY(160) 

SION  RF0A( 201 ) , RF0B( 201 ) , RFOC( 201 ) , S( 201 ) 

SION  OMAA(201),OMBB(201),OMCC(201),TZ(160) 

SION  W(  201  )  ,  WW(  201  )  ,  TV.'(  1 60)  ,  TX  (  1 60 )  ,  TY  ( 1 60) 

SION  X1(160),Y1(160),Z1(160),X2(160),Y2(160) 

SION  X3(160),Y3(160),Z3(160),X4(160),Y4(160) 

SION  X0(160),Y0(160),Z0(160),Z4(160) 

SION  FIL0U1(606),  FIL0U2 ( 606 ) , FIL0U3( 730) 

ed  Mxy ( 1 , 2 , 3 ) , IDE=1  from  calcu.  curve;  =2  saved 
egree(l),  90  degree(3);  IDG=2  no  optimization; 
,  Dxy(2);  ID1I  =  4  taking  half  of  the  180  as  90; 
11,*)  DD.DT.II.ML.IDP 

0M,M,C0C,P,NL 

11,12,13, 14, GI.SIKC 

DRTA , RFA , COE , COEF , DRTA 1 

IDD, IDE, IDF, IDG, IDH, IDO 


11 
11 
11 
11 


READ(11 
PRINT  * 
PRINT  * 
PRINT  * 
PRINT  * 
PRINT  * 
PRINT  * 
PRINT  * 


*) 
*) 
*) 
*) 


*)  KC.P.ATIO.KD 
DD,DT,H,ML,IDP 
OM,M,COC,P,NL 
I1,I2,I3,I4,GI,SINC 

DRTA , RFA , COE , COEF , DRTA 1 
IDD .IDE, IDF, IDG, IDH, IDO 
RXY1,RXY2,RY1,RY2,RZ1,RZ2 
KC, RATIO, KB 


N3  =  N 

;;2=3 
if  ( 

READ 
IF  ( 
DO  7 
RFOA 
RFOB 
RFOC 


L-l 
*N3 
IDE. 
(10, 
IDH. 
0  1  = 
(D  = 
(D  = 
(D  = 


EQ.2)  THEN 

*)  (IDUM.FILINP(L) 

EQ.4)  THEN 

1.N3 

FILINP(I)/2 

FILINP(I+N3)/2 

FILINP(I+2*N3)/2 


L=1,M2) 
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70  CONTINUE 
ELSE 

DO  71  1=1, N3 
RFOA(I)=FILINP(I) 
RF0B(I)=FILINP(I+N3) 
RF0C(I)=FILINP(I+2*N3) 

71  CONTINUE 
ENDIF 

ENDIF 

:a  =  (NL-i)/2+i 

N2=N1-1 

N4=Nl-2 

RN1=N1 

Iil  =  (HL+l)/2 

RM1=M1 

IF  (IDO.EQ.l)  THEN 

DO  72  1=1,11 
DESFIL(I)=H 

72  CONTINUE 

IF  (IDD.EQ.l)  THEN 

IF  (IDF.E0.3)  THEN 

DESFIL(I1+1)=0.5*H 

DESFIL(I1+I2+2)=0.5*H 

ELSE 

DESFIL(I1+1==0.0 

DESFIL(Il+I2+2)=0.0 

ENDIF 

DO  73  1=1,12 

DESFIL(I+I1+1)=-GI 

73  CONTINUE 

DO  74  1=1,13 
DESFIL(I+I1+I2+2)=H 

74  CONTINUE 
ENDIF 

IF  (IDD.EQ.2)  THEN 

IF  (IDF.EQ.3)  THEN 

DESFIL(Il+l)=3*H/4 

DESFIL(Il  +  2)=II/2 

DESFIL(Il+3)=H/4 

DESFIL(Il+3+I2+l)=H/4 

DESFIL(Il+3+I2+2)=H/2 

DESFIL( I 1+3+12+3 )=3*H/4 

ELSE 

DESFIL(Il+l)=H/2 

DESFIL(Il+2)=0.0 

DESFIL(Il+3)=-GI/2 

DESFIL( I 1+3+12+1 )=-GI/2 

DESFIL(Il+3+I2+2)=0.0 

DESFIL(Il+3+I2+3)=H/2 

ENDIF 

DO  75  1=1,12 
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DESFIL(Il+3+I)=-GI 
7  5  CONTINUE 

DO  76  1=1,13 

DESFIL(Il  +  I2  +  6  +  I)=II 
7  6  CONTINUE 

ENDIF 

IF  (IDD.EQ.3)  THEN 

DESFIL(Il+l)=3*GI/4 

DESFIL(Il+2)=GI/2 

DESFIL(Il+3)=GI/4 

DESFIL(Il+4)=0.0 

DESFIL(Il+5)=-GI/4 

DESFIL(Il+6)=-GI/2 

DESFIL(Il+7)=-3*GI/4 

DESFIL(Il+7+I2+l)=-3*GI/4 

DESFIL(Il+7+I2+2)=-GI/2 

DESFIL( I 1+7+12+3 )=-GI/4 

DESFIL(Il+7+I2+4)=0.0 

DESFIL(Il+7+I2+5)=GI/4 

DESFIL(Il+7+I2+6)=GI/2 

DESFIL( I 1+7+12+7 )=3*GI/4 

DO  77  1=1,12 

DESFL(Il+7+I)=-GI 
7  7  CONTINUE 

DO  78  1=1,13 

DESFIL( I 1+7+1 2+7+1 )=GI 

78  CONTINUE 
ENDIF 
ENDIF 

C 

IF  (ID0.EQ.2)  THEN 
DO  79  1=1,11 
DESFIL(I+3*HL)=0.0 

79  CONTINUE 

C  desired  profile  is  Mxy 
IF  (IDD.EQ.l)  THEN 
DESFIL(I1+1+3*ML)=0.5*H 
DESFIL(Il+I2+2+3*ML)=0.5*H 
DO  80  1=1,12 
DESFIL(I1+1+I+3*ML)=H 

80  CONTINUE 

DO  81  1=1,13 
DESFIL(Il+I2+2+I+3*ML)=0.0 

81  CONTINUE 

ELSEIF  (IDD.EQ.2)  THEN 

DESFIL(1+I1+3*ML)=0.25*H 

DESFIL(2+I1+3*KL)=0.5*H 

DESFIL(3+I1+3*ML)=0.75*H 

DESFIL(l+Il+I2+3+3*ML)=0.7  5*H 

DESFIL(2+Il+I2+3+3*HL)=0.5*H 

DESFIL(3+Il+I2+3+3*ML)=0.25*H 

DO  32  1=1,12 

DESFIL(3+I1+I+3*KL)=H 

82  CONTINUE 
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DO  83  1=1,13 
DESFIL(I+6+Il+I2+3*ML)=0.0 

83  CONTINUE 
ENDIF 

i 

ENDIF 

DO  84  I  =  l,3*iIL 

DESFIL(I)=0.0 

84  CONTINUE 

URITE(1S,1200)  (L,PESFIL(L),L=1,4*KL) 
1200  F0RiIAT(I5,E14.6) 
DO  35  1=1, ML 
DXY(I)=DESFIL(I+3*ML) 
8  5  CONTINUE 

i 

IK  =  1 

112  BETA=0.0 
K  =  l 

DO  86  N=1,NL 
S(N)=0.0 
S1(N)=0.0 
86  CONTINUE 

111  DO  100  L-l.ML 

RL=L 

IF  (L.LE.M1)  THEN 

D0M-(RL-RM1)*DD 

ELSEIF  (L.GT.M1)  THEM 

D0M=(RL-RM1)*DD 

ENDIF 

X(1)=0.0 

Y(1)=0.0 

Z(1)=H 

DO  101  N=1,N3 

rn=n 

IF  (K.EQ.l)  THEN 

IF  (IDE.EQ.l)  THEN 

IF  (M.EQ.l)  THEN 

0MA=0M 

OME=OM 

OMC=OM 

elseif  (m.eq.2)  then 
exa=((rn-rn1)*dt)**2/(2*d**2) 

o;;a=om*exp(-exa) 

EXB=((RN-RNl)*DT+DT/2)**2/(2*D**2) 

0MB=0M*EXP(-EXB) 

EXC=((RN-PN1)*DT+PT)**2/(2*D**2) 

OMC=OM*EXP(-EXC) 

ELSEIF  (M.EQ.3)  THEN 

IF  (N.EQ.N1)  THEN 

OHA=OM 

ELSE 

exa=p*(rn-rn1)*dt 
oi;a=oh*sin(exa)/(exa) 
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ENDIF 

EXB=P*((RN-  :Nl)*DT+DT/2) 

OMB-OM*SI  (  :XB)/(EXB) 

if  (N.EQ.N2)  the;; 
o:;c=om 

ELSE 

exc=p*((en-rni)*dt+dt) 
omc=om*sin(exc)/(exc) 

ENDIF 

ENDIF 

ELSEIF  (IDE.EQ.2)  THEN 

OMA=RFOA(N) 

OMB=RFOB(N) 

OMC=RFOC(N) 

ENDIF 

OMAA(N)-OMA 

OMBB(N)=OMB 

OMCC(N)-OMC 

ELSEIF  (K.NE.l)  THEN 

oi;a=omaa(m) 
omb=ombb(n) 
omc-omcc(n) 

ENDIF 

IF  (K.EQ.l)  THEN 

RF1A(N)=0HA 

RF1B(N)=0MB 

RF1C(N)=0MC 

ENDIF 

Al=DT*DOM*Y(N) 
Bl=DT*(-DOM*X(N)-OMA*Z(N)) 
Cl=DT*OMA*Y(N) 
A2=DT*DON*(Y(N)+Bl/2) 

B2=DT*(-DOM*(X(N)+Al/2)-OKB*(Z(N)+Cl/2)) 
C2=DT*OMB* ( Y ( N ) +B1 / 2 ) 
A3=DT*DOH*(Y(N)+B2/2) 
B3=DT*(-DOM*(X(N)+A2/2)-OMB*(Z(N)+C2/2)) 

C3=DT*0;iE*(Y(N)  +  B2/2) 
A4=DT*DOM*(Y(N)+B3) 

B4=DT*(-DOM*(X(N)+A3)-OMC*(Z(N)+C3)) 
C4=DT*0MC*(Y(N)+B3) 
X(N+l)=X(N)+(Al+2*A2+2*A3+A4)/6 
Y(N+l)-Y(N)+(Bl+2*B2+2*B3+B4)/6 
Z(N+l)-Z(N)+(Cl+2*C2+2*C3+C4)/6 
101  CONTINUE 

IF  (K.EQ.l)  THEN 

IF  (L.E0.M1)  THEN 

PRINT  *,  Z(1),Z(N1),Z(NL) 

ENDIF 

ENDIF 

ZO(L)=Z(NL) 

XO(L)=X(NL) 
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YO(L)=Y(NL) 
100  CONTINUE 

IF  (IDF.EQ.3)  THEN 

DO  102  L=1,ML 

RL  =  L 

D0M=-(RL-P.M1)*DD*0.5 

X(1)=X0(L) 

Y(D  =  YO(L) 

Z(1)=Z0(L) 

W(1)=0.0 

DO  103  N-1.N3 

OIIA  =  0.0 

Oi!B=0.0 

0MC=0.0 

A1=DT*D0M*Y(N) 

B 1=DT* ( -DOM*X ( N ) -OMA*Z ( N ) ) 

C1=DT*0MA*Y(N) 

D1=D1*(X(N)**2+Y(N)**2-DXY(L)**2)**2 

A2=DT*D0M*(Y(N)+Bl/2) 

B2=DT*(-D0M*(X(N)+Al/2)-0MB*(Z(N)+Cl/2)) 

C2=DT*OMB*(Y(N)+Bl/2) 

DD= ( Y ( N )+B 1*0.5) **2-DXY ( L) **2 

D2=DT*((X(N)+A1*0.5)**2+DD)**2 

A3=DT*DOM*(Y(N)+B2/2) 

B3=DT*(-D0M*(X(N)+A2/2)-0MB*(Z(N)+C2/2)) 

C3=DT*0MB* ( Y( N )+B2/2 ) 

DD=(Y(N)+B2*0.5)**2-DXY(L)**2 

D3-DT*((X(N)+A2*0.5)**2+DD)**2 

A4=DT*D0M*(Y(N)+B3) 

B4=DT*(-D0M*(X(N)+A3)-0MC*(Z(N)+C3)) 

C4=DT*0MC*(Y(N)+B3) 

DD=(Y(N)+B3)**2-DXY(L)**2 

D4=DT*((X(N)+A3)**2+DD)**2 

X(N+l)«X(N)+(Al+2*A2+2*A3+A4)/6 

Y(N+l)=Y(N)+(Bl+2*B2+2*B3+B4)/6 

Z(N+l)=Z(N)+(Cl+2*C2+2*C3+C4)/6 

W(N+l)=W(N)+(Dl+2*D2+2*D3+D4)/6 
103  CONTINUE 

IF  (K.EO.l)  THEN 

Z1(L)=Z(NL) 

X1(L)=X(NL) 

Y1(L)=Y(NL) 

ENDIF 

TX(L)=X(NL) 

TY(L)=Y(NL) 

TW(L)=W(NL) 
102  CONTINUE 

ENDIF 

J0=0.0 
JA=0.0 
DO  21  L=1,KL 


c 
c 
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J=RXY2*(TX(L)**2+TY(L)**2-DXY(L)**2)**2+TW(L) 

JO=JO+J 
J=RXY2*(TX(L)**2+TY(L)**2-DXY(L)**2)**2 

JA=JA+J 
21  CONTINUE 

IF  (K.EQ.l)  THEN 
PRINT  *,  JO.JA 
ENDIF 
no  optimization,  goto  to  the  end 
IF  (IDG.EQ.2)  THEN 
GOTO  402 
ENDIF 

DO  200  L=1,ML 

RL=L 

D0M=(RL-RH1)*DD 

XA(1)=0.0 

YA(1)=0.0 

ZA(1)=H 

DO  201  N=1,N3 

RK-N 

IF  (K.EQ.l)  THEN 

IF  (IDE.EQ.l)  THEN 

IF  (M.EQ.l)  THEN 

0MA=0H 

OMB-OH 

0MC=0H 

ELSEIF  (M.EQ.2)  THEN 

EXA=((RN-RN1)*DT)**2/(2*D**2) 

OMA=OM*EXP(-EXA) 

EXB-((RN-RNl)*DT+DT/2)**2/(2*D**2) 

0MB=0h*EXP(-EXB) 

exc=((rn-rn1)*dt+dt)**2/(2*d**2) 

o;jc=om*exp(-exc) 
elseif  (m.eq.3)  then 
if  (n.eq.n1)  then 
oka=om 

ELSE 

EXA=P*(RN-RN1)*DT 
OMA=OM*SIN(EXA)/(EXA) 
ENDIF 

EXB=P*((RN-RNl)*DT+DT/2) 

OMB=OM*SIN(EXB)/(EXB) 

IF  (N.EQ.N2)  THEN 

OMC=OM 

ELSE 

EXC=P*((RN-RN1)*DTDDT) 

OMC=OM*SIN(EXC)/(EXC) 

ENDIF 

ENDIF 

ELSEIF  (IDE.EQ.2)  THEN 
0MA=RFOA(N) 
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OMB 

o;:c 

END 

ELS 
OMA 
OMB 
OMC 
END 

integra 
Al  = 
Bl  = 
Cl  = 
A2  = 
B2  = 
C2  = 
A3  = 
33  = 
C3  = 
A  4  = 
B4= 
C4  = 
XA( 
YA( 
ZA( 
CON 
RL= 
DOM 
XB( 
YB( 
ZB( 
W(l 

integra 
DO 
OMA 


!01 


=RFOB(N) 

-RFOC(N) 

IF 

EIF  (K.NE. 
=OMAA(N) 
=OMBB(N) 
=OMCC(N) 
IF 

te  Eq.[66] 
DT*D0i:*YA( 
DT*(-DOM*X 
DT*OMA*YA( 
DT*DOM*(YA 
DT*(-DOM*( 
DT*OMB*(YA 
DT*DOM*(YA 
DT*(-DOM*( 
DT*OMB*(YA 
DT*DOM*(YA 
DT*(-DOM*( 
DT*OMC*(YA 
N+1)=XA(N) 
N+1)=YA(N) 
N+1)=ZA(N) 
TINUE 
L 

=-(RL-RMl) 
1)=X0(L) 
1)«Y0(L) 
1)=Z0(L) 
)  =  0.0 

te  Eq.[72] 
211  N-1.N3 
=  0.0 
OMB=0.0 
OMC=0.0 
A1=DT*D0M*YB( 
B1=DT*(-D0M*X 
C1=DT*0MA*YB( 
B1=DT*(XB(N)* 
A2»DT*D0K*(YB 
32=DT*(-D0M*( 

C2=dt*o;:b*(yb 

DD=(YB(N)+B1* 
D2=DT*((XB(N) 
3=DT*D0M*(YB( 
B3=DT*(-D0M*( 
C3=DT*0MB*(YB 
DD=(YB(N)+B2* 
DB=D?*((XE(N) 
A4=DT*D0K*(YB 
B4=DT*(-D0M*( 
C4=DT*0HC*(YE 
I)D=(YB(N)+B3) 


1)  THEN 


of  Chap 
N) 

A(N)-OKA 
N) 

(TO+B1/2 
XA(N)+A1 
(IJ)  +  Bl/2 
(N)+B2/2 
XA(N)+A2 
(N)+B2/2 
(N)+B3) 
XA(N)+A3 
(N)+B3) 
+(A1+2*A 
+(B1+2*B 
+(C1+2*C 


*DB*0.5 


ter  V 
*ZA(N)) 

) 
/2)-0MB*(ZA(N)+Cl/2)) 

) 
) 
/2)-OMB*(ZA(N)+C2/2)) 

) 
)-0MC*(ZA(N)+C3)) 

2+2*A3+A4)/6 

2+2*B3+B4)/6 
2+2*C3+C4)/6 


of  Chapter  V 


N) 

B(N)-OMA 

N) 

*2+YB(M) 
(N)+Bl/2 
XBfN)+Al 
(IO  +  B1/2 
0.5)**2- 
+A1*0.5) 
N)+B2/2) 
XE(K)+A2 

(N)+B2/2 
0.5)**2- 

+A2*0.5) 

(N)+B3) 

XB(N)+A3 

(N)+B3) 

**2-DXY( 


*ZB(N)) 

**2-DXY(L)**2)**2 

) 
/2)-0MB*(ZB(N)+Cl/2)) 

) 
DXY(L)**2 

**2+DD)**2 

/2)-OMB*(ZB(N)+C2/2)) 

) 

DXY(L)**2 

**2+DD)**2 

)-0MC*(ZE(N)+C3)) 

L)**2 


115 


D4=DT*((B(N)+A3)**2+DD)**2 
XE(N+l)»XB(N)+(Al+2*A2+2*A3+A4)/6 
YB(N+l)=»YB(N)  +  (Bl+2*E2+2*B3+B4)/6 
ZB(N+l)=ZB(N)+(Cl+2*C2+2*C3+C4)/6 
W(N+l)=!v'(N)  +  (Dl+2*D2+2*D3+D4)/6 
211  CONTINUE 
C  integrate  Eq.[84]  of  Chapter  V 
X  ( 1 )  =  1 .  0 
Y(1)=0.0 
Z(1)=0.0 
U(1)=0.0 
BO  202  N=1,N3 
O!!A=0.0 
OMB-0.0 

o;ic=o.o 

XD=(XE(N)+XB(N+l))/2 
YD=(YB(N)+YB(N+l))/2 
Al=DT*DOH*Y(N) 
B1=DT*(-D0M*X(N)-0MA*Z(N)) 
C1=DT*0MA*Y(N) 

DA=2*(XB(N)**2+YB( N ) **2-DXY ( L ) **2 ) *2*XB ( N ) *X ( H ) 
DB=2*(XB(N)**2+YB(N)**2-BXY(L)**2)*2*YB(N)*Y(N) 
D1=DT*(DA+DB) 
A2-DT*DOM*(Y(N)+Bl/2) 

B2=DT*(-D0M*(X(N)+Al/2)-0MB*(Z(N)+Cl/2)) 
C2-DT*OHB*(Y(N)+Bl/2) 

DA=2*(XD**2+YD**2-DXY(L)**2)*2*XD*(X(N)+A1*0.5) 
DB=2*(XD**2+YD**2-DXY(L)**2)*'2*YD*(Y(N)+B1*0.5) 
D2=DT*(DA+DB) 
A3=DT*DOM*(Y(N)+B2/2) 

B3=DT* ( -DOM* (X(N)+A2/2 ) -OKB* ( Z( N )+C2/2 ) ) 
C3=DT*OMB*(Y(N)+E2/2) 

DA=2*(XD**2+YD**2-DXY(L)**2)*2*XD*(X(N)+A2*0.5) 
DB=2* ( XD**2+YD**2-DXY( L ) **2 )*2*YD* ( Y ( N )+B2*0 . 5 ) 
D3  =  DT*(1)A  +  DB) 
A4-DT*D0H*(Y(N)+B3) 

B4-DT*  ( -DOM*  ( X  (  N  )  +A3  )  -CMC*  (  7,  (  N  )  +C3 ) ) 
C4=DT*0MC*(Y(N)+B3) 
AA=2*XB( N+l ) * ( X (N)  + A3 ) 
BB=2* YB ( N+l ) * ( Y( N ) +B3 ) 

DA=2*(XB(N+1)**2+YB(N+1)**2-DXY(L)**2)*AA 
DB=2* (XB(N+1 ) **2+YB (N+l )**2-DXY (L)**2) *BB 
D4oDT*(DA+DB) 

X(N+l)=X(N)+(Al+2*A2+2*A3+A4)/6 
Y(K+l)=Y(N)+(Bl+2*B2+2*B3+B4)/6 
Z(N  +  l)=Z(lI)  +  (Cl  +  2*C2+2*C3  +  C4)/6 
W(N+l)=W(N)  +  (Dl+2*D"2+2*D3+D4)/6 
202  CONTINUE 

Z3(L)=Z(NL) 

X3(L)=X(NE) 

Y3(L)=Y(NL) 

W3(L)«W(NL) 

DZ=X3(L)*(TX(L)**2+TY(L)**2-DXY(L)**2)*2*TX(L) 

EZ=Y3(L)*(TX(L)**2+TY(L)**2-DXY(L)**2)*2*TY(L) 
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XX 
C  i  n  t  e  g  r 
X( 
Y( 
Z( 
W( 

do 

OH 
OH 
O'l 
XD 
YD 
Al 
Bl 
CI 
DA 
DB 
Dl 
A2 
E2 
C2 
DA 
DB 
D2 
A3 
B3 
C3 
DA 
DD 
D3 
A  A 
FA 
C4 
A  A 
BB 
DA 
DB 
DA 
X( 


(1)-DZ+EZ+W 

ate  Eq.[84] 


3(L) 
of  Chapter  V 


212 


Y( 
Z( 
W( 

CO 
Z3 
X3 
Y3 


DZ 


1)=0.0 
1)-1.0 
1)=0.0 
l)=0.0 

212  M=1,K3 
A  =  0.0 
B=0.0 
C=0.0 

=(XB(N)+XB 
=(YB(N)+YB 
=DT*D0M*Y( 
=DT*(-DOH* 
=DT*0MA*Y( 
=2*(XB(N)* 
=2*(XE(N)* 
=DT*(DA+DB 
«DT*DOM*(Y 
-DT*(-DOM* 
=DT*OMB*(Y 
=2*(XD**2+ 
=2*(XD**2+ 
=DT*(DA+DB 
=DT*DOM*(Y 
=DT*(-DOM* 
=DT*OMB*(Y 
=2*(XD**2+ 
«2*(XD**2+ 
=DT*(DA+DE 
=DT*DOM*(Y 
-DT*(-DOM* 
=DT*OMC*(Y 
-2*XB(K+1) 
=2*YB(N+1) 
=2*(XB(N+1 
=  2*(XB(NT+1 
=DT*(DA+DB 
M+1)=X(N)+ 
N+1)-Y(N)+ 
N+1)=Z(N)+ 
N+1)-W(N)+ 
NTINUE 
(L)-Z(NL) 
(L)-X(NL) 
(L)-Y(NL) 
(L)-W(NL) 
=X3(L)*(TX( 


N+l))/2 

N+l))/2 

) 
(N)-OMA*Z(N)) 

) 

2+YB ( N ) **2-DX Y ( L) **2 )*2*XB( H ) *X ( N  ) 

2+YB(M ) **2-DXY (L)**2) *2*YB(N)*Y(N ) 

N)+Bl/2) 

X(N)+Al/2)-0MB*(Z(X)+Cl/2)) 

N)+Bl/2) 

D**2-DXY( L ) ** 2 ) *2*XD* (X (N )+A 1 *0 . 5 ) 
D**2-DXY( L) **2 ) *2*YB* ( Y(N )+Bl*0 . 5 ) 

N)+52/2) 

X ( N ) +A2/ 2 ) -OMB* (Z(N )+C2/2 ) ) 

N)+B2/2) 

D**2-DXY(L)**2)*2*XD*(X(N)+A2*0.5) 

D**2-DXY ( L ) * *  2 ) *  2 * YD* ( Y ( N ) +B2*0 . 5 ) 

N)+B3) 

X ( N ) +A3 ) -OMC* (Z (N )+C3) ) 

N)+B3) 

(X(N)+A3) 

(Y(N)+B3) 

**2+YB( N+ 1 ) **2-DXY ( L ) **2 ) * A A 

**2+YB(N+l)**2-DXY(L)**2)*BB 

A1+2*A2+2*A3+A4)/G 
Bl+2*B2+2*B3+B4)/6 
Cl+2*C2+2*C3+C4)/6 
Dl+2*D2+2*D3+D4)/6 


YY 
C  integr 


L)**2TY(L)**2-DXY(L)**2)*2*TX(L) 
EZ=Y3(L)*(TX(L)**2+TY(L)**2-DXY(L)**2)*2*TY(L) 

3(L) 
of  Chapter  V 


(1)=DZ+EZ+W; 

ate  Eq.fSA] 
X(l)-0.6 
Y(1)=0.0 

Z(l)=1.0 
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222 


C 
C 


W(l)- 

DO  22 

o:;a=o 

OMB=0 

o;;c=o 

XD  =  (X 
YD=(Y 
A1  =  DT 
B1=DT 
C1  =  DT 
DA  =  2* 
DB=2* 
D1=DT 
A2=DT 
B2=DT 
C2  =  DT 
DA=2* 
DB=2Ji 
D2=DT 
A3=DT 
B3  =  DT 
C3=DT 
DA=2* 
DB=2* 
D3  =  DT 
A4  =  DT 
B4=DT 
C4  =  DT 
AA=2* 
BB=2* 
DA-2* 
DB=2* 
D4=D* 
X(N+1 
Y(N+1 
Z(N+1 

w(n+i 

CONTI 
Z3(L) 
X3(L) 
Y3(L) 
W3(L) 
DZ=X3 
EZ=Y3 
ZZ(1) 


0.0 

2    K-l.N 

.0 

.0 

.0 

B(N)+XB 

B(N)+YB 

*DOM*Y( 

*(-DOM* 

*OMA*Y( 

(XB(N)* 

(XB(K)* 

*(DA+DE 

*DOM*(Y 

*(-DOM* 

*OMB*(Y 

(XD**2+ 

(XD**2+ 

*(DA+DB 

*DOM*(Y 

*(-dom* 

*OMB*(Y 
(XD**2+ 
(XD**2+ 
*(DA+DB 
*DOM*(Y 
*(-DOM* 
*OMC*(Y 
XB(N+1) 
YB(N+1) 

(X3(: +i 

(XB(N+1 
(DA+DB) 
)=X(N)+ 
)=Y(N)+ 
)=Z(N)+ 
)=W(N)+ 

MI1T7 

11    U  Jj 

-Z(NL) 
=X(NL) 

=Y(NL) 

=W(NL) 

(L)*(TX 

(L)*(TX 

=DZ+EZ+ 


N+l))/2 
N+l))/2 

) 
(K)-OMA*Z(N)) 

) 

!2+YB(N)**2-DXY(L): 
!2+YB(N)**2-DXY(L)! 


*2)*2*XB(N)*X(N) 
*2)*2*YB(N)*Y(N) 


;;)+bi/2) 

X(N)+Al/2)-OMB*(Z(N)+Cl/2) ) 
N)  +  1:1/2) 

D**2-DXY(L)**2)*2 
D**2-DXY(L)**2)*2 


*XD*(X(N)+A1*0.5) 
*YD*(Y(N)+B1*0.5) 


)+E2/2) 
X ( N )  +  A  2 / 2 ) -OMB* ( Z ( N ) +C2 / 2 ) ) 
N)+B2/2) 

D**2-DXY( L ) **2 ) *  2 
D**2-DXY(L)**2)*2 


*XD*(X(N)+A2*0.5) 
*YD*(Y(N)+B2*0.5) 


,:)+B3) 

X(N)+A3)-0KC*(Z(N)+C3)) 

N)+B3) 

(X(N)+A3) 

(Y(N)+B3) 

**2+YE(N+l)**2-DX 

**2+YB(N+l)**2-DX 


(A1+2*A2+2*A3+A4)/ 
(B1+2*B2+2*B3+B4)/ 
(C1+2*C2+2*C3+C4)/ 

(D1+2*D2+2*D3+D4)/ 


Y(L)**2)*AA 
Y(L)**2)*BB 

6 
fc> 
6 
6 


(L)**2+TY(L)**2-DXY(L)**2)*2*TX(L) 

( L ) **2+TY ( L ) **2-DXY ( L ) **2 )*2*TY( L ) 
W3(L) 


1)0M=-2*D0M 


DO 

203    N  = 

-1, 

K3 

RN  = 

-.11 

DD1 

:=-dt 

IF 

(K.EQ 

.1) 

TI 

IEN 

IF 

(IDE.] 

::o. 

15 

THEN 

IF 

(H.EO 

■  i) 

TI 

IEN 
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OMA-OM 

OMB-OM 

OKC-OM 

ELSEIF  (M.EQ.2)  THEN 

EXA=((RN-RN1)*DT)**2/(2*D**2) 

0MA=0M*EXP(-EXA) 

EXB=((RN-RNl)*DT+DT/2)**2/(2*D**2) 

OMB=OM*EXP(-EXB) 

EXC=((RN-RN1)*DT+DT)**2/(2*D**2) 

OMC=OM*EXP(-EXC) 

ELSEIF  (M.EQ.3)  THEN 

IF  (N.EQ.N1)  THEN 

OMA-OM 

ELSE 

EXA=P*(RN-RN1)*DT 

OMA=OM*SIN(EXA)/(EXA) 

END  IF 

EXE=P*((RN-RNl)*DT+DT/2) 

OMB-OM*SIN ( EXB ) / ( EXB ) 

IF  (N.EQ.N2)  THEN 

OMC-OM 

ELSE 

EXC=P*((RN-RN1)*DT+DT) 

OMC-OM* S I N ( EXC ) / ( EXC ) 

ENDIF 

ENDIF 

ELSEIF  (IDE.EQ.2)  THEN 
OMA-RFOA(N) 
OMB-RFOB(N) 
OMC-RFOC(N) 
ENDIF 

ELSEIF  (K.NE.l)  THEN 
OMA-OMAA(N) 
OMB-OMBB(N) 
OMC-OMCC(N) 
ENDIF 
integrate  Eq.[92]  backwards  in  time 
Al=DDT*DOM*YY(N) 
Bl»DDT*(-DOM*XX(N)-OMA*ZZ(N)) 
C1=DBT*0MA*YY(N) 
A2=DDT*DOM*(YY(N)+Bl/2) 

B2-DDT*(-D0M*(XX(N)+Al/2)-OHB*(ZZ(N)+Cl/2)) 
C2-DDT*OIIB*  (  YY  (  K  )  +B 1  /  2  ) 
A3=DUT*DOM*(YY(N)+B2/2) 

B3=DDT*(-DOM*(XX(N)+A2/2)-OMB*(ZZ(N)+C2/2)) 
C3-DDT*0MB*(YY(N)+B2/2) 
A4=DDT*D0M*(YY(N)+B3) 

B4=DDT*(-D0I1*(XX(N)  +  A3)-0MC*(ZZ(N)+C3)) 
C4=DDT*OMC*(YY(N)+B3) 
XX(N+l)=XX(N)+(Al+2*A2+2*A3+A4)/6 
YY(N+l)=YY(N)+(Bl+2*B2+2*B3+B4)/6 
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ZZ(N+l)«ZZ(N)+(Cl+2*C2+2*C3+C4)/6 

203  CONTINUE 

DO  2  3  K-l.NL 

G( M) =-ZA (N ) *YY ( ML+ 1-N ) + Y A ( N ) *ZZ( NL+ 1-N) 

23  CONTINUE 

IF  (L.EQ.l)  THEN 
DO  2  4  N-l.NL 
GG(N)-0.0 

24  CONTINUE 
ENDIF 

DO  25  N-l.KL 
GG(N)=GG(N)+G(N) 

25  CONTINUE 

200  CONTINUE 
SG1-0.0 
DO  31  1=1, N2 
SG1=SG1+(C0E*GG(2*I-1))**2 

31  CONTINUE 
SG2=0.0 

DO  32  1=1, N4 
SG2=SG2+(C0E*GG(2*I))**2 

32  CONTINUE 

SG= ( C0E*GG( 1 ) ) **2+( COE*GG ( NL ) ) **2 
D0T2=(SG+4*SGl+2*SG2)*DT/3 
IF  (K.NE.l)  THEN 
BETA=D0T2/D0T1 
ENDIF 
D0T1=D0T2 
DO  33  N=1,NL 
S ( N ) =-GG ( N ) +BETA*S ( N ) 
3  3  CONTINUE 

PRINT  *,  GG(1),GG(N1),GC(NL),K 
PRINT  *,  S(1),S(N1),S(NL),BETA 

HUA=G. 613034 

RFA1=RFA*HA 

RFA2=RFA 

L1=RFA*HUA**2 

L2=RFA*HUA 

L3=0.0 

L4=RFA 

RLN-L1-L3 

RLM=L4-L2 

IF  (K.NE.l)  THEN 

IF  (K.GT.IDP)  THEM 

RFA=RFA*C0C 

ENDIF 

RFA1=C0EF*RFA*HUA 

RFA2=C0EF*RFA 

L1=C0EF*RFA*HUA**2 

L2=C0EF*RFA*HUA 

L3=0.0 
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333 


L4=C0EF*RFA 

RLN=L1- 

■L3 

RLM=L4- 

-L2 

ENDIF 

ID  =  0 

1  =  0 

DO  300 

L-l.ML 

RL=L 

D0M=(RL-RM1)*DD 

X(1)=0. 

0 

Y(1)=0. 

0 

Z(l)-H 

DO  301 

N=l ,N3 

RN=N 

IF  (K.EQ.l)  THEN 

IF  (IDE.EQ.l)  THEN 

IF  (M.EQ.l)  THEN 

0MA=0M+RFA1*S(N) 

0;iB=0M+RFAl*(S(N)  +  S(N  +  l))/2 

0MC=0M+RFA1*S(N+1) 

ELSEIF  (M.EQ.2)  THEN 
EXA=((RN-RN1)*DT)**2/(2*D**2) 

ONA=OH*EXP(-EXA)+RFA1*S(N) 
EXB=((RN-RNl)*DT+DT/2)**2/(2*D**2) 

0MB=0M*EXP(-EXB)+RFAl*(S(N)+S(N+l))/2 
EXC=((RN-RN1)*DT+DT)**2/(2*D**2) 
0iiC  =  0H*EXP(-EXC)+RFAl*S(N  +  l) 

ELSEIF  (M.EQ.3)  THEN 

IF  (N.EQ.N1)  THEN 

0MA=0M+RFA1*S(N) 

ELSE 

EXA=P*(RN-RN1)*DT 

0MA=0M*SIN(EXA)/(EXA)+RFA1*S(N) 

ENDIF 

EXB=P*((RN-RNl)*DT+DT/2) 

0ME=0K*SIN(EXB)/(EXB)+RFAl*(S(N)+S(N+l))/2 

IF  (N.EQ.N2)  THEN 

0MC=0M 

ELSE 

EXC=P*((RN-RN1)*DT+DT) 

0MC=0M*SIN(EXC)/(EXC)+RFA1*S(N+1) 

ENDIF 

ENDIF 

ELSEIF  (IDE.EQ.2)  THEN 

0KA=RF0A(N)+RFA1*S(N) 

OMB=RF0B(N)+RFAl*(S(N)+S(++l))/2 

0MC=RF0C(N)+RFA1*S(N+1) 

ENDIF 
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ELSEIF  (K.NE. 
0MA=0ilAA(N)+R 
OMB=OI;BB(N)  +  R 
OMC-OKCC(N)+R 
ENDIF 
integrate  Eq . [ 66 ] 

o;:aa(n)=oma 
ombb(n)»omb 

ohcc(n)=omc 

RF2A(N)=0MA 

RF2B(N)=0MB 

RF2C(N)=0MC 

A1»DT*D0M*Y(N 

B1=DT*(-D0M*X 

C1=DT*0HA*Y(N 

A2=DT*D0M*(Y( 

B2*DT*(-D0M*( 

C2=DT*0MB*(Y( 

A3=DT*D0M*(Y( 

B3-DT*(-D0M*( 

C3=DT*0KB*(Y( 

AA=DT*DOM*(Y( 

B4*=DT*(-D0M*( 

C4=DT*0HC*(Y( 

X(N+l)=X(N)+( 

Y(N+l)=Y(N)+( 

Z(N+l)=Z(N)+( 

301  CONTINUE 

ZO(L)=Z(NL) 
XO(L)=X(NL) 
YO(L)=Y(NL) 

300  CONTINUE 


1)  THEN 

FA1*S(N) 

FAl*(S(N)+S(N+l))/2 

FA1*S(N+1) 

of  Chapter  V 


) 
(N)-0MA*Z(N)) 

) 


N) 
X( 
N) 
N) 
X( 
N) 
N) 
X( 
N) 
Al 
El 
CI 


+  B1 
N)  + 

+  B1 
+  B2 

R)  + 

+  32 
+  B3 
N)  + 
+  B3 
+  2* 
+2* 
+  2* 


/2) 

Al/2)-0MB*(Z(N)+Cl/2)) 

/2) 

/2) 

A2/2)-OMB*(Z(N)+C2/2)) 

/2) 

) 

A3)-0MC*(Z(N)+C3)) 

) 

A2+2*A3+A4)/6 
B2+2*B3+B4)/6 
C2+2*C3+C4)/6 


IF  (IDF.EQ.3)  THEN 
DO  302  L-l.ML 
RL=L 

D0M=-(RL-RM1)*DD*0.5 
X(1)=X0(L) 
Y(1)=Y0(L) 
Z(1)=Z0(L) 
W(1)=0.0 
integrate  E q . [ 7 2 ]  of  Chapter  V 
DO  303  N=1,N3 
01  i  A  =  0.0 
0MB=0.0 

onc=o.o 

Al=DT*DOM*Y(N 
Bl=DT*(-DOM*X 
Cl-DT*OMA*Y(N 
D1=D1*RXY1*(X 
A2«DT*D0M*(Y( 
B2=DT*(-D0M*( 
C2=DT*OMB*(Y( 
DD=(Y(N)+B1*0 


) 
(N)-0MA*Z(N)) 

) 

(N)**2+Y(N)**2-DXY(L)**2)**2 

N)+Bl/2) 

X(N)  +  Al/2)-0JIB*(Z(N)+Cl/2)) 

IO  +  B1/2) 

,5)**2-DXY(L)**2 
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D2=DT*RXY1*((X(N)+A1*0.5)**2+DD)**2 

A3=DT*D0M*(Y(N)+B2/2) 

B3=DT*(-DOM*(X(N)+A2/2)-OMB*(Z(N)+C2/2)) 

C3=DT*0MB*(Y(N)+32/2) 

DD= ( Y( N )+B2*0 . 5 ) #*2-DXY(L)**2 

H3»DT*RXY1*((X(N)+A2*0.5)**2+DD)**2 

A4=DT*D0M*(Y(N)+B3) 

B4=BT* ( -DOM* ( X ( N ) +A3 ) -OMC* ( Z ( N ) +C3 ) ) 

C4=DT*0MC*(Y(N)+B3) 

DD=(Y(N)+B3)**2-DXY(L)**2 

D4=DT*RXY1 * ( ( X ( N ) +A3 ) **  2  +  DD ) **2 

X(N+l)«X(N)+(Al+2*A2+2*A3+A4)/6 

Y(N+l)=Y(N)+(Bl+2*B2+2*B3+B4)/6 

Z(N+l)=Z(N)+(Cl+2*C2+2*C3+C4)/6 

V(N+l)=W(N)+(Dl+2*D2+2*D3+D4)/6 
303  CONTINUE 

Z2(L)=Z(NL) 

X2(L)=X(IIL) 

Y2(L)=Y(NL) 

TW(L)-W(NL) 

TX(L)=X(NL) 

TY(L)-Y(NL) 
302  CONTINUE 

ENDIF 

JJ=0.0 

JB=0.0 

DO  41  L=1,HL 
J*=RXY2*(TX(L)**2+TY(L)**2-DXY(L)**2)**2+TW(L) 

JJ=JJ+J 

J  =  RXY  2* ( TX ( L )**2+TY ( L ) **2-DXY (L) **2 )**2 
JB=JB+J 
41  CONTINUE 
1  =  1  +  1 

IF  (I.EQ.l)  THEN 
J2=JJ 

RFA1-RFA1*HUA 
GOTO  333 

ELSEIF  (I.EQ.2)  THEN 
J1=JJ 
ENDIF 

IF  (ID.EQ.l)  THEN 
J1=JJ 

ELSEIF  (ID.EQ.2)  THEN 
J2=JJ 
ENDIF 

IF  (K.GT.l)  THEN 

IF  (Jl.GT.JC)  THEN 

IDO-1 

GOTO  400 

ELSEIF  (JB.GT.JD)  THEN 

ID0=1 

GOTO  400 
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EX  DIP 
ENDIF 

IF  (K.EQ.I4)  THEN 

ID0=1 

GOTO  400 

ENDIF 

IF  (J2.LT.J1)  TEEN 

RFA2=RFA2-RLN 

ID=2 

L3=L1 

L1=L2 

J1=J2 

L2=L3+HUA*RFA2 

RLM-L4-L2 

RLN=L1-L3 

EFA1=L2 

GOTO  333 

ELSEIF  (J2.GT.J1+DRTA)  THEN 

RFA2=RFA2-RLH 

ID=1 

L4=L2 

L2=L1 

J2=J1 

L1=L4-HUA*RFA2 

RLM=L4-L2 

RLN=L1-L3 

RFA1=L1 

GOTO  333 

ELSE 

JC=J1 

JD=JB 

PRINT  *,  0,J1,JA,JB 

ENDIF 

DO  43  N=1,LL 
S1(N)=S1(N)+RFA1*S(N) 
4  3  CONTINUE 
K  =  K+1 

IF  (J2.LT.J1+DRTA1)  THEN 
ID0=1 
GOTO  400 
ELSE 
ID0=2 
GOTO  111 
ENDIF 
400  CONTINUE 

IF  (IDE.EQ.2)  THEN 
output  as  input,  run  again 
IF  (IK.LT.KD)  THEN 
DO  401  1=1, N3 
RF0A(I)=RF2A(I) 
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RF0B(I)=RF2B(I) 

RF0C(I)=RF2C(I) 

401  CONTINUE 
IK=IK+1 
PRINT  *,IK 
GOTO  112 
ENDIF 
ENDIF 

i 

402  DO  44  KK=1,N3 
FILGU1(KK)=RF1A(KK) 
FIL0U1 (KK+N3 ) =RF1E ( KK ) 
FIL0U1(KK+2*N3)=RF1C(KK) 

44  CONTINUE 

DO  45  KK-1.N3 

FILOU2(KK)=RF2A(KK) 

FILOU2(KK+N3)=RF2B(KK) 

FIL0U2(KK+2*N3)=RF2C(KK) 
4  5  CONTINUE 

DO  46  KK=1,ML 

FIL0U3(KK)=Z1(KK) 

FIL0U3 ( KK+ML ) =X 1 ( KK ) 

FIL0U3(KK+2*ML)=Y1(KK) 

SQ=X1(KK)**2+Y1(KK)**2 

FIL0U3(KK+3*ML)=SQRT(SQ) 

46  CONTINUE 

DO  47  KK11.ML 

FIL0U4(KX)=Z2(KK) 

FIL0U4(KK+ML)=X2(KK) 

FILOU4(KK+2*ML)=Y2(KK) 

SQ=X2(NK)**2  +  Y2(KK)':*2 

FIL0U4(KK+3*HL)=SQRT(SQ) 

47  CONTINUE 
M3-4*ML 

DO  48  KK=1,M2 
FIL0U5(KK)=0.0 

48  CONTINUE 

DO  49  KK=1,K3 
FILOU6(KK)=0.0 

49  CONTINUE 

WRITE(22,1100)  (L,FIL0U1(L),L=1,M2) 
WRITE  (2 1,1100)  (L,FILOU2(L),L=l,!I2) 
WRITE( 14 , 1 100)  ( L , FIL0U3 ( L) , L=l , M3 ) 
WRITE  (15, 11 00)  (L,FIL0U4(L),L=  1,113) 
V;RITE(16,1100)  (L,FIL0U5(L),L=1,M2) 
WRITE(17,1100)  (L,FILOU6(L),L=l,M3) 

C  suppress  the  floating  point  underflow 
CALL  ERRSET(208, 256,-1, 1) 

1100  FORMAT (1 5, E44. 6) 
STOP 
END 


APPENDIX  I 
JOB  PROGRAM  FOR  APPENDICES  B  AND  H 

THE  JOB  FOR  APPENDIX  E: 

//MAOl  JOB  (2003,5001, 100, 7, 999), MAO, CLASS=1 

//  EXEC  FORTVCLE 

/♦INCLUDE  SELE51  FORTRAN 

//G0.FT18F001  DD  SYSOUT=C 

//G0.FT10F001  DD  * 

/♦INCLUDE  FILINP  DATA 

//G0.FT10F001  DD  * 

/♦INCLUDE  INPI 

//G0.FT22F001 

//GO.FT13F001 

//G0.FT14F001 

//G0.FT15F001 

//G0.FT16F001 

//G0.FT17F001 

//GO.SYSIN  DD 

THE  JOB  FOR  APPENDIX  H: 

//MA02  JOB  (2003,5001, 100, 7, 999), MAO, CLASS-1 

//  EXEC  FORTVCLE 

/♦INCLUDE  SELE52  FORTRAN 

//G0.FT18F001  DD  SYSOUT=C 

//G0.FT10F001  DD  ♦ 

/♦INCLUDE  FILINP  DATA 

//G0.FT10F001  DD  ♦ 

/♦INCLUDE  INPI 

//G0.FT22F001 

//G0.FT13F001 

//GO.FT14F001 

//G0.FT15F001 

//G0.FT16F001 

//G0.FT17F001 

//GO.SYSIN  DD 


TF 

DATA 

DD 

SYSOUT=J 

DD 

SYSOUT-E 

DD 

SYSOUT=F 

DD 

SYSOUT=G 

DD 

SYSOUT-K 

DD 

SYSOUT-I 

* 

TF 

DD 

DATA 
SYSOUT=J 

DD 

SYSOUT=E 

DD 

SYSOUT=F 

DD 

SY30UT=G 

DD 

SYSOUT=H 

DD 

SYSOUT=I 
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APPENDIX  J 
PARAMETER  VALUES  OF  THE  PROGRAM  IN  APPENDIX  H  TO  OBTAIN 
TEE  OPTIMAL  RF  PULSES  OF  FIG.  20,  FIG.  25  AND  FIG.  28 

a)  The  parameter  values  for  the  optical  rf  pulse  of  Fig.  20 


1570.0, 

1.11111 1E-05 , 

IOC, 

159, 

2 

-9065.0, 

0.7, 

19300.0, 

91 

66, 

21, 

66, 

5, 

o, 

1 

50.0, 

0.1E-03, 

0.5E-05, 

0.1E-03, 

10.0 

2, 

1, 

3, 

0, 

5, 

2 

1, 

1, 

1, 

1, 

1, 

1 

20, 

0.7, 

15 

b )  The  parameter  values  for  the  optimal  rf  pulse  of  Fig.  25 


1570.0, 

1.1] 

11111E-05, 

100, 

159, 

1 

-9065.0, 

0.7, 

19300.0, 

91 

67, 

23 

67, 

5, 

o, 

1 

50.0, 

0.1E-03, 

0.5E-05, 

0.1E-03, 

10.0 

1, 

2, 

3, 

o, 

5, 

2 

1, 

1, 

1, 

1, 

1, 

1 

20, 

0.7, 

12 

c)  The  parameter  values  for  the  optimal  rf  pulse  of  Fig.  28 


1570.0, 

1, 

.  11111 1E-05, 

100, 

159, 

10 

-15396.0, 

3, 

0.7, 

31975.0, 

91 

58, 

37, 

53, 

5, 

o, 

1 

50.0, 

0.1E-03, 

0.5E-05, 

0.1E-03, 

10.0 

2, 

1, 

3, 

0, 

5, 

2 

1, 

1, 

1, 

1, 

1, 

1 

20, 

0.7, 

15 
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APPENDIX  K 
FOURIER  COEFFICIENTS  OF  THE  FOURIER  SERIES  REPRESENTATION 
OF  THE  OPTIMAL  RF  PULSE  IN  FIG.  20 

Sine 

O.OOOOOOE+OO 
-0.201899E+03 

0.507699E+03 
-0.337984E+03 
-0.383514E+00 
-0.579742E+02 
-0.342352E+G2 
-0.407180E+02 
-0.337707E+02 
-0.316671E+02 
-0.279499E+02 
-0.251290E+02 
-0.233887E+02 
-0.208346E+02 
-0.192330E+02 
-0.182061E+02 
-0.155616E+02 
-0.1586S3E+02 
-0.134981E+02 
-0.130264E+02 
-0.125601E+02 
-0.106989E+02 
-0.109925E+02 
-0.951433E+01 
-0.912414E+01 
-0.867219E+01 
-0.771478E+01 
-0.740869E+01 
-0.689631E+01 
-0.605017E+01 
-0.626775E+01 
-0.489746E+01 
-0.523200E+01 
-0.427739E+01 
-0.39S430E+01 
-0.390304E+01 
-0.303956E+01 
-0.304002E+01 
-0.249087E+01 
-0.203356E+01 
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Cosine 

1 

-0. 

, 169372E+04 

2 

0. 

.288795E+04 

3 

-0, 

.361280E+04 

4 

0. 

.159014E+04 

5 

0. 

,141494E+01 

6 

0. 

.159420E+03 

7 

0. 

.769519E+02 

8 

0. 

, 765642E+02 

9 

0. 

, 540168E+02 

10 

0. 

■435509E+02 

11 

0. 

.332784E+02 

12 

0. 

, 259571E+02 

13 

0. 

.210232E+02 

14 

0. 

.162115E+02 

15 

c. 

•129513E+02 

16 

0. 

.104777E+02 

17 

0. 

.755129E+01 

18 

0. 

.633576E+01 

19 

0. 

.434571E+01 

20 

0. 

.325267E+01 

21 

0. 

.216887E+01 

22 

0, 

•110160E+01 

23 

0. 

•369705E+00 

24 

-0, 

.350803E+00 

25 

-0, 

.974517E+00 

26 

-0. 

.156207E+01 

27 

-0. 

.193531E+01 

28 

-0. 

.242C71E+01 

29 

-0, 

.278906E+01 

30 

-0. 

•298010E+01 

31 

-0. 

.364296E+01 

32 

-0. 

.332352E+01 

33 

-0, 

.415392E+01 

34 

-0, 

.383264E+01 

35 

-0, 

.415360E+01 

36 

-0, 

.464790E+01 

37 

-0, 

.422144E+01 

38 

-0, 

.486577E+01 

39 

-0, 

.469502E+01 

40 

-0, 

.457746E+01 

APPENDIX  L 
FOURIER  COEFFICIENTS  OF  THE  FOURIER  SERIES  REPRESENTATION 
OF  THE  OPTIMAL  RF  PULSE  IN  FIG.  28 

Sine 

O.OOOOOOE+OO 
-0.199309E+03 

0.447691E+G3 
-0.616315E+03 

0.103004E+04 
-0.620669E+04 

0.413190E+02 
-0.114660E+03 
-0.552481E+02 
-0.567666E+02 
-0.403002E+02 
-0.386490E+02 
-0.310853E+02 
-0.291230E+02 
-O.257630E+02 
-0.230600F+02 
-0.209952E+02 
-0.190140E+02 
-0.177844E+02 
-0. 162040E+02 
-0.151512E+02 
-0.139326E+02 
-0.123477E+02 
-0.121618E+02 
-0.111424E+02 
-0.10496&E+02 
-0.955069E+01 
-0.899339E+01 
-0.828694E+01 
-0.757017E+01 
-0.727276E+01 
-0.634151E+01 
-0.602479E+01 
-0.542441E+01 
-0.500013E+01 
-0.437506E+01 
-0.333065E+01 
-0.357354E+01 
-0.292428E+01 
-0.264080E+01 
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Cosine 

1 

-0. 

158867E+04 

2 

0. 

285033E+04 

3 

-0, 

318531E+04 

4 

0. 

290026E+04 

5 

-0. 

359238E+04 

6 

0. 

170538E+04 

7 

-0. 

934551E+02 

8 

0. 

214829E+03 

9 

0. 

830845E+02 

10 

0. 

778256E+02 

11 

0. 

47S703E+02 

12 

0. 

398490E+02 

13 

0. 

27891 1E+02 

14 

0. 

226311E+02 

15 

0. 

173047E+02 

16 

0. 

132494E+02 

17 

0. 

101745E+02 

18 

0. 

762904E+01 

10 

0. 

571372E+01 

2  0 

0. 

402426E+01 

21 

0. 

262465E+01 

22 

0. 

144144E+01 

23 

0. 

424740E+00 

24 

-0. 

445334E+00 

25 

-0. 

120196E+01 

26 

-0. 

190271E+01 

27 

-0. 

240136E+01 

28 

-0 

295330E+01 

29 

-0. 

336101E+01 

30 

-0. 

370644E+01 

31 

-0  . 

422037E+01 

32 

-0. 

430422E+01 

33 

-0. 

473654E+01 

34 

-0. 

434606E+01 

35 

-0. 

518969E+01 

36 

-0 

522512E+01 

37 

-0 

531313E+01 

38 

-0 

573393E+01 

39 

-0 

549935E+01 

40 

-0 

596654E+01 
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